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The evolution of galaxy cluster X-ray scaling relations 
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ABSTRACT 

We use numerical simulations to investigate, for the first time, the joint effect of feedback 
from supernovae (SNe) and active galactic nuclei (AGN) on the evolution of galaxy cluster 
X-ray scaling relations. Our simulations are drawn from the Millennium Gas Project and 
are some of the largest hydrodynamical iV-body simulations ever carried out. Feedback is 
implemented using a hybrid scheme, where the energy input into intracluster gas by SNe and 
AGN is taken from a semi-analytic model of galaxy formation. This ensures that the source 
of feedback is a population of galaxies that closely resembles that found in the real Universe. 
We show that our feedback model is capable of reproducing observed local X-ray scaling 
laws, at least for non-cool core clusters, but that almost identical results can be obtained 
with a simplistic preheating model. However, we demonstrate that the two models predict 
opposing evolutionary behaviour. We have examined whether the evolution predicted by our 
feedback model is compatible with observations of high-redshift clusters. Broadly speaking, 
we find that the data seems to favour the feedback model for z < 0.5, and the preheating 
model at higher redshift. However, a statistically meaningful comparison with observations is 
impossible, because the large samples of high-redshift clusters currently available are prone 
to strong selection biases. As the observational picture becomes clearer in the near future, it 
should be possible to place tight constraints on the evolution of the scaling laws, providing us 
with an invaluable probe of the physical processes operating in galaxy clusters. 

Key words: hydrodynamics - methods: N-body simulations - galaxies: clusters: general - 
galaxies: cooling flows - X-rays: galaxies: clusters. 



1 INTRODUCTION 

Galaxy cluster surveys are a potentially powerful means of plac- 
ing tight constraints on key cosmological parameters, independent 
of other methods such as the measurement of cosmic microwave 
background anisotropics. This is primarily because the mass func- 
tion of galaxy clusters is highly sensitive to different choices of 
model parameters. It is therefore essential to determine how the 
mass function varies with redshift if we are to exploit clusters as 
a cosmological probe. This is not trivial since the total masses of 
galaxy clusters must first be inferred from their observable proper- 
ties. 

Mass estimates can be derived from X-ray observations of the 
hot, diffuse intracluster medium (ICM) by assuming that the gas 
is in hydrostat ic equilibrium within the cluster gravitational poten- 
tial well (e.g. ISarazirJ 1988). However, this requires the accurate 
determination of gas density and temperature profiles out to large 
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radii. Although this has become common practice for low-redshift 
clusters with the advent of Chandra and XMM-Newton, measuring 
high-quality profiles of distant clusters requires long observation 
times, so hydrostatic mass estimates have only been made for a 
small number of high-redshift clusters. Furthermore, the assump- 
tion of hydrostatic equilibrium only applies to dynamically relaxed 
systems, so this technique cannot be applied to unbiased cluster 
samples. 

In cases where a hydrostatic estimate is not possible, clus- 
ter masses can be inferred from the relationships that exist be- 
tween X-ray observables, like luminosity, Lx, temperature, T, and 
the total mass, M. These scaling relations are predicted by the 
simple self-similar model of cluster formation, where the ICM is 
heated solely by gravitational processes, such as adiab atic compres - 
sion and shocks induced by supersonic accretion jKaiseilll986t) . 
However, the observed Lx-M relat ion is steeper than expected 
from gravitationa l heating alone (e.g-j Reiprich & Bohringerl l2002l : 
IChen etalJ l200l IVikhlinin et al] 120091: Pratt et alj |200aC as is 
the Lx-T relation (e.g.lMarkevitcri[l998l;lAjnaud & Evrardlll9"99l 
IWuet alJI 19991 : lEttori et alJl2002l ; |Pratt et alj|2009h - This departure 
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from self-similarity i s due to an excess of entropvj in cluster cores 
dPonman et alj|l999l ; Illovd-Davies et alj|200 Oj; iFinoguenov et al.l 
2002). This extra entropy prevents gas from being compressed to 
high densities, reducing its X-ray emissivity compared to the self- 
similar prediction. The effect will be more pronounced in galaxy 
groups since they have shallower potential wells, leading to a steep- 
ening of the Lx-M relation (and Lx-T relation) as desired. The 
main physical processes thought to be responsible for boosting the 
entropy of the ICM are heating from astrophysical sources, such 
as supernovae (SNe) and active galactic nuclei (A GN), and the re- 
moval of low-entropy gas via radiative cooling (see Voit 2005] for a 
review). 

Unfortunately, scaling relations are not a perfect means of 
converting from X-ray observables to mass because of their in- 
trinsic scatter. The dominant source of scatter about relations in- 
volving Z/x is the dense, hi ghly X-ray luminous central regions o f 
cool core (CC) clusters (e.g. lo Hara et alj|2006l ; lchen et al.1l2007l) . 
One way of reducing this scatter is simply to exclud e the core 
region from the measurements (e.g. Mark evitchiri998l) . Another 
source of scatter is cluster mergers, which can cause clusters to 
shift (approxi mately) along the Lx-T relation, but awa y from the 
M-T relation dRowlev et al.l2004 l Kravtsov et alj2006l) . However, 
iKravtsov et alT ?2006 ) recently demonstrated that the quantity Yx, 
defined as the product of the gas mass Af gas and the X-ray temper- 
ature, is extremely tightly correlated with total mass and is insensi- 
tive to cluster mergers. This result has been c onfirmed by indepen- 
dent numerical simulations ( Poo le et al .12007]) and several observa- 



tional studies dArnaud et al 



2007; Maughan 



2007; Vikhlinin etai] 



2009), suggesting that reliable cluster mass estimates can indeed be 
derived from simple observables. 

Given that scaling relations enable us to estimate the masses 
of clusters with poor quality X-ray data, it is clearly vital to en- 
sure that they are well calibrated over a wide redshift range, so we 
can improve cosmological constraints derived from the mass func- 
tion. Furthermore, measuring the evolution of X-ray scaling laws 
should reveal information on the nature of the non-gravitational 
heat ing and cooling process es responsible for shaping galaxy clus- 
ters dMuanwong et alj |2006). 

At low redshift, scaling relations are reasonably well cali- 
brated, at least for relaxed clusters, since reliable hydrostatic mass 
estimates are readily available. Measuring these relations at high- 
redshift is considerably more challenging, with fewer results avail- 
able in the literature. As will become clear from the following brief 
revi ew, the observationa l picture is far from settled at present. 

lEttori et alj d2004bl) used a sample of 28 galaxy clusters with 
0.4 < z < 1.3 to measure the M-T, Lx-M and Lx-T rela- 
tions. Upon comparing their high-redshift M-T relation with local 
counterparts, they found that the normalisation evolved in a manner 
well described by the self-similar model. By contrast, the normal- 
isation of both the high-redshift Lx-M and Lx-T relations was 
lower than expected from self-similar scaling argume nts. Negative 
evolut ion^ of the Lx-T relation was also reported bv lHilton et al.l 
d2007l) , based on XMM-Newton observatio ns of a single cluster a t 
z = 1.457. Another study was performed by Maughan et al.l 120061) 



1 We define entropy as K = k^T / n2 _1 , where fee is Boltzmann's con- 
stant, n c is the electron number density and 7 = 5/3 is the ratio of specific 
heats for a monoatomic ideal gas. 

2 Throughout this work, we consider the evolution of scaling relations rel- 
ative to that predicted by the self-similar model. Given a particular scaling 
law, we say there is negative (positive) evolution if the normalisation at 
high-redshift is lower (higher) than anticipated from self-similar scaling. 



who used a sample of 11 clusters in the redshift range 0.6 < z < 1. 
They concluded that there is no evidence for any evolution of the 
M-T, Lx-M and Lx-T relations beyond that anticipated from 
self-simi lar theory. Similar resu lts f or the Lx-T 1 relatio n were also 
found bv lVikhlinin et all d2002l) and lLumb etafl (2004). However, 
such analyses often rely on simplifying assumptions, like isother- 
mal temperature profiles and /3-model profiles for the gas den- 
sity, which can lead to large b iases in the mass determination (e.g. 
iMarkevitch & VikhlinirJI 19971) . 

A different approach was adopted by Morandietal. (2007), 
whose sample consisted of 24 clusters in the redshift range 0.14 < 
z < 0.82. For each object in their sample, they used the measured 
gas density profile and the equation of hydrostatic equilibrium to 
determine the best-fitting temperature profile by varying the free 
parameters in the assumed model for the dark matter density pro- 
file. In conflict with most of the results discussed above, they re- 
ported slight negative evolution of the M-T and Lx-M relations, 
and positive evolution of the Lx-T relation. 

More recently, small samples of distant clusters (~ 10 ob- 
jects) with hydrostatic mas s estim ates have been used to measure 
the Yx-M (Maughanl2007h. M-T dKotov & Vikhlinirj|2005l.l200r3) 
and Lx-T dKotov & Vikhlininl 20051) relations at high-redshift (z ~ 
0.7). The evolution of the normalisation of the Yx-M and M-T re- 
lations was again found to be consistent with the self-similar model. 
On the other hand, the Lx- T relation exhibited p ositive evolution, 
consistent with the results of lMorandi et aljj2 007). but not the neg- 
ative e volution measured by lEttori et al.l d2004bl) and iHilton et al.l 
d2007l) . Beyond z ~ 0.7, there are very few clusters with a hydro- 
static mass estimate; a notable exc eption is the z = 1.05 object dis- 
covered by Mau ghan et alj d2008bl) . who demonstrated that the po- 
sition of this cluster on the high-redshift Yx-M, M-T and Lx-M 
relations can be explained by self-similar scaling arguments. 

Although several observational studies indicate that the self- 
similar model provides an adequate description of the evolution of 
X-ray scaling relations, IVoitl d2005l) suggests that self-similar evo- 
lution cannot persist to arbitrarily high redshift, because of the in- 
creasing importance of radiative cooling and feedback from galaxy 
formation. An exa mple of an observational result that supports this 
argument is that of iBranchesi et alj d2007l) . Using a sample of 39 
distant clusters (0.25 < 2 < 1.3), they found that the evolution 
of the Lx-T relation is comparable with the self-similar prediction 
for < z < 0.3, but negative at higher redshifts. 

One reason for the lack of concordance between different ob- 
servational studies is inconsistent cluster selection. With current 
heterogeneous archival cluster samples, it is very difficult to ac- 
count for selection biases, which may vary from sample to sample 
and can mimic evolution. An attempt to quantify the impact of se- 
lection effects on th e evolution of scaling relations was made by 
IPacaud et alj d2007h , who focused on the Lx-T relation. From an 
analysis of their raw data, they found evidence for non-monotonic 
evolution, at odds with the self-similar prediction. However, after 
taking selection effects into account, they found that the evolution 
of the Lx-T relation was consistent with self-similar scaling, al- 
though this was a tentative result because of poor statistics. Never- 
theless, the salient point to take from their work is that modelling 
the full source-selection process can drastically alter the measured 
evolution, so we must seriously question any attempt to assess the 
evolution of X-ra y scaling laws that do es not attempt to do this. In 
more recent work. lMantz et alj J2009bl lah have developed a method 
for rigorously accounting for selection effects in order to constrain 
the scaling relations at low and high redshifts, including their evolu- 
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tion. They concluded that the T-M and Lx-M relations both show 
no departures from self-similar evolution up to z ss 0.5. 

In the near future, the XMM Cluster Survey (XCS; 
iRomer et alj200l|) will provide a complete sample of clusters span- 
ning a broad redshift range, < z < 1.5, all coming from the same 
survey and selected with well-defined criteria. Since the survey se- 
lection function will be known, the evolution of X-ray scaling re- 
lations will be measured with unprecedented accuracy out to high- 
redshift. Cluster surveys based on the Sunyaev-Zel'dovich (SZ) ef- 
fect, such as those bein g conducted with the South Pole Telescope 
(Vandeiiind e et alj|2010h and the Atacama Cosmology Telescope 
( IHincks et alj|2009h . are also expected to produce large catalogues 
of clusters covering a wide redshift range with a well-defined mass 
threshold. Like XCS, these samples will be ideal for studying the 
evolution of cluster scaling relations. 

There are also several other reasons why results from differ- 
ent studies may appear to be contradictory. For example, the mea- 
sured evolution can be affected by: poor statistics due to small 
sample sizes; the local scaling relation chosen to compare high- 
redshift data with; different definitions of evolution (i.e. whether 
the expected self-similar behaviour is scaled out first or not); dif- 
ferent choices of scale radius used to define a cluster (i.e. redshift- 
dependent or not); and whether the core region is excised or not 
(and the size of the core region). In short, extreme care is required 
when comparing observational results with each other, and with 
theoretical results from numerical simulations. 

Theoretical studies of cluster formation have mainly con- 
centrated on explaining the lack of self-similarity inherent in 
low-redshift scaling relations, rather than the evolution of scal- 
ing laws. To this end, numerous mechanisms for raising the 
entropy of the intracluster gas have been implemented in hy- 
drody namica l TV-body simulations, such as radiative coolin. 
(e.g. iBrvarJ l200d: iPearce et alj l200d; iMuanwong et alj 1200 J 



|2002| ; iDave et al.1 120021: IVoit et alj I200I; IWu & Xud |2002|). pre 
heating (e.g. iBialek et alj 1200 J iBrighenti & MathewsTl200ll ; 



Muanwong et alj|2002l : iBorgani et alj|2002l ; iTornatore et alj|2003l ; 
Borga ni et aL 20051). star formation and a ssociated feedback from 



Valdarnini 2003; Kav 2004|; iKav et"al] |2004 2007; 



2007; Borgani et al. 2004, 2005; R omeo et alj |2006) 



and black hole grow t h with associated feedback from AGN ( e.g. 
iPuchwein et al.ll2008l : lMcCarthv et al.ll2010l : lFabian et alfeoiol) . 

Simulations employing such models have indeed proved capa- 
ble of successfully reproducing the observed scaling relations for 
local clusters, at least on average, although often at the expense of 
excessive star formation. However, we would not expect the ther- 
mal history of the ICM to be the same in each case, so we should see 
di fferences in evolutionary behaviour. This was first demonstrated 
bv lMuanwong et all j2006h who traced the evolution of the cluster 
population to z — 1.5 using three different models: a cooling-only 
model, a model with preheating and cooling, and a stellar feedback 
model that self-consistently heats cold gas in proportion to the local 
star formation rate. While all three schemes produce indistinguish- 
able Lx-T relations at z = 0, they predict strongly positive, mildly 
positive and mildly negative evolution of the Lx-T relation, re- 
spectively. Therefore, given observational data of sufficient quality, 
it should be possible to rule out unsuitable theoretical models of 
non-gravitational heating. This emphasises the importance of accu- 
rately measuring cluster scaling relations over a range of redshifts. 

Other attempts to study the evol ution of X-ray scali ng laws 
with numerical simulations are scarce. lEttori et ail (2004a) used a 
simulation including radiative cooling, star formation and super- 
nova feedback to follow the evolution of the cluster population be- 



tween z = 1 and z = 0.5. They measured a small, but significant, 
negative evolution in the normalisation of the Lx-M and Lx-T re- 
lations, and a marginally negative evolution of the M-T relation. 
These results were found t o be consistent w ith the observational 
data of lEttori et all <2004bh . IKav et al.1 <2007h used a similar simu- 
lation, but with a different prescription for cooling and stellar feed- 
back, to investigate the evolution of the M-T and Lx-T relations 
in the redshift interval < z < 1. The results the y obtained are 
qualitatively the same as those o flEttorietalJj2004iJ) . However, we 
note that a rigorous comparison bet ween the results o f lEttori et al.1 
d2004iJ) , lMuanwong et alJ d~2006) and Kav et aljj2007t) is not possi- 
ble owing to differences in their data analysis procedures (different 
definitions of cluster scale radii, etc.) 

In this paper, we reconsider the evolution of cluster scaling re- 
lations from a theoretical perspective. Our basic goal is to investi- 
gate how scaling laws evolve when feedback from AGN is included 
in numerical simulations, in addition to stellar feedback. As far as 
we are aware, this is the first time this has been attempted. 

The simulation we use is drawn from the Millennium Gas 
Project, a series of hydrodynamical A^-body simulations with 
the same volume (500 3 /i -3 Mpc 3 ) and the same initial per- 
turbation amplitudes and phases as the Millennium Simulation 
dSpringel et alj|2005l) . Feedback from g alaxies is incorporated i n 
our simulation via the hybrid approach of lShort & Thomas! J2009h . 
where the energy imparted to the ICM by SNe and AGN is com- 
puted from a semi-analytic model (SAM) of galaxy formation. 
Since SAMs are tuned to reproduce the properties of observed 
galaxies, the source of feedback in our simulation is guaranteed to 
be a realistic galaxy population. By contrast, fully self-consistent 
simulations with radiative cooling, star formation and supernova 
feedb ack typically predict that too muc h gas cools and forms stars 
(e.g.Hprgani et al. 2004; Kav et al. 200%. We assess how feedback 
from galaxy formation affects the evolution of scaling relations by 
comparing our results with those obtained from two other Millen- 
nium Gas simulations. The first of these is a control model, where 
the gas is heated by gravitational processes only. The second in- 
cludes additional radiative cooling and uniform preheating at high 
redshift as a simple model of non-gravitational heating from astro- 
physical sources. 

The large volume of the Millennium Gas simulations enables 
us to resolve statistically significant numbers of clusters at all red- 
shifts relevant to cluster formation. In fact, our cluster samples are 
some of the largest ever extracted from numerical simulations. For 
example, in the preheating plus cooling run we have ove r 20 times 
more o bjects with T > 2 keV at z — and z = 1 than lKav et al.1 
d2007l) . Furthermore, we can capture the formation of the richest 
systems, allowing us to probe the same dynamic range as the ob- 
servations. The Millennium Gas simulations thus provide an ideal 
tool for studying the evolution of the cluster population. 

The layout of this paper is as follows. In Section[2]we briefly 
review the self-similar model of cluster formation. The three Mil- 
lennium Gas simulations are discussed in Section [3] along with 
a description of our method for generating cluster catalogues and 
profiles. In Section|4]we compare results for our three models, start- 
ing with a discussion of cluster profiles and X-ray scaling relations 
at z = 0, then investigating how profiles and scaling laws evolve 
from z = 1.5 to z = 0. Wherever possible we attempt to place 
observational constraints on our models. We summarise our results 
and conclude in Section[5] 
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2 CLUSTER SCALING THEORY 

The simplest model of galaxy cluster formation is that they form 
via the gravitational collapse of the most overdense regions in the 
dark matter distribution, and the cluster baryons are heated only 
by gravitational processes (compression and shock heating) during 
the collapse. Since non-linear gravitational processes do not intro- 
duce any characteristic scale, we would then expect clusters to be 
self-similar, i.e. scaled versions of each other. With the additional 
assumptions that clusters are spherically symmetric systems and 
that the intracluster gas is in hydrostatic equilibrium with the un- 
derlying dark matter potential, it is straightforward to deri ve sim- 
ple se lf-similar scaling relations between cluster properties jKaised 
U9M). 

Defining r a as the radius of a spherical volume within which 
the mean matter density is A times the critical density at redshift z, 
we find that the total enclosed mass, 

MaE 

scales with gas temper- 
ature, Ta, as 

E( z y 2/3 T A oc M A /S . (1) 

Under the further assumption that the X-ray emission of the ICM 
is primarily thermal Bremsstrahlung radiation (which is valid for 
T > 2 keV), the luminosity, Lx,a, within r A is given by 

EOO^Lx.AOcTi. (2) 

The scaling between X-ray luminosity and total mass follows upon 
combining equations QJ and (|2}: 

E{z)~ 7/3 Lx,a oc M A /3 . (3) 

The quantity Yx = M gaa T has recently attracted much atten- 
tion since it has been shown to be a l ow-scatter mass proxy , regard- 
less of cluster dynamical state (e.g. iKravtsov et al]|2006l) . This is 
primarily because Yx approximates the total thermal e nergy of the 
ICM, which is not stron gly affected by cluster me rgers dPoole et all 
120071) . unlike L x or T dRicker & Sarazinll200lh . The self-similar 
scaling of Yx within ta with total mass is 

£(*r 2/3 Yx,A oc M A /3 . (4) 

The density contrast A governs the scale radius within which 
one measures the mass of a cluster. The most common choice is 
to set A = 500, since rsoo is the effective limiting radius for re- 
liable observations from Chandra and XMM-Newton. Throughout 
this paper we will adopt A = 500, independent of redshift. 

Note that there are other ways in which th e scale radiu s of a 
cluster can be defined. In the original model o f Kaiserl fl98^> . clus- 
ters were assumed to be self-similar with respect to the mean mat- 
ter density of the Universe, rather than the critical density. In this 
case, the E(z) factors in equations (fTJ-lO are replaced by (1 + z). 
Another possibility is to assume a redshift-dependent density con- 
trast A(z) that is proportional to A v i r , z , where A v j r ( z)p CI (z) is 
the mean density within the cluster virial radius (e .g. lEttori et all 
l2004b! : iMaughan et alj|2006l ; iBranchesi et alj|2007h . The value of 
Avir(2) is taken from the analytical solution for the collapse of a 
spherical top-hat perturbation, under the assumption th at the clus- 
ter has just virialised at that redshift (e.g. lPeebles 1980). With this 
choice of density contrast, the E(z) factors in equations {TJ — ([4j» 
become E(z) A(z) 1 ^ 2 instead. Therefore, predictions for the evo- 
lution of the cluster population clearly depend on how the scale 

3 M A = 47rr A Ap cr (^)/3, where p CI (z) = 3H$E(z) 2 /SwG is the 
critical density and E(z) 2 = f2 m o(l + z) 3 + ^A,o m a spatially-flat 
ACDM cosmological model. 



radius is defined, so caution must be exercised when comparing the 
results of different observational and theoretical studies. 

Observations of clusters in the local Universe have established 
that scaling relations between cluster properties do indeed exist, 
but their form is found to be different to the self-similar predic- 
tions. This is because non-gravitational effects, such as radiative 
cooling and feedback from galaxy formation, modify the entropy 
structure of the ICM. Lower-mass systems are affected more than 
massive objects, breaking the self-similarity of the scaling laws. 
Departures from self-similarity thus provide us with a probe of 
the non-gravitational processes operating in clusters. Although the 
self-similar model cannot explain the observed form of scaling re- 
lations, there is some observational evidence that it is capable of 
describing their evolution. 



3 SIMULATIONS: THE MILLENNIUM GAS PROJECT 

In this work we use simulations taken from the Millennium Gas 
Project, a progra mme to add gas to t he dark matter-only Millen- 
nium Simulation dSpringel e t al. 2005). We present a new addition 
to the Millennium Gas suite in which feedback is directly tied to 
galaxy formation, enabling us to investigate how energy input from 
both SNe and AGN affects the evolution of cluster scaling rela- 
tions. Hereafter, we refer to this simulation as the FO run. The 
feedba ck model adopted is the hybrid scheme of IShort & T homas 
(2009), where a SAM is used to calculate the energy transferred 
to the intracluster gas by SNe and AGN. An immediate benefit 
of this approach is that feedback is guaranteed to originate from a 
galaxy population whose observational properties agree well with 
those of real galaxies. This is generally not the case in fully self- 
consistent hydrodynamical simulations that include radiative cool- 
ing and stellar feedback because too much gas cools out of the hot 
phase, leading t o excessive star formation (e.g. lBorgani et alj2004l : 
iKav et al.ll2007l) . It is widely thought that additional heating from 
AGN is the natural so lutio n to this ove rcooling problem. Indeed, 
iMcCarthv etai] feo id) and lFabian etai] feOlfJ) have demonstrated 
that including AGN feedback in hydrodynamical simulations can 
successfully balance radiative cooling in galaxy groups. However, 
the stellar fraction is still fo und to be ~ 2 — 3 times larger than 
observed in massive clusters (Fabjan et al. 2010h . 

By coupling a SAM to a hydrodynamical simulation, 
IShort & Tho mas (2009) were able to reproduce the observed mean 
Lx-T relation for groups and poor clusters, provided there was 
a large energy input into the ICM from AGN over the entire for- 
mation history of haloes. The AGN heating is more efficient at 
driving X-ray emitting gas from the central regions of lower-mass 
haloes, reducing their luminosity and steepening the Lx-T rela- 
tion as desired. Their model was also able to account for some of 
the increased scatter about the mean relation seen for temperatures 
T < 3 keV, attributable to the varie d merger histories of gro ups. 

A limitation of the method of IShort & Thomas! d2009l) is that 
radiative cooling is not incorporated in their hydrodynamical sim- 
ulations. Instead, their model relies on the simplistic treatment of 
the distribution and cooling of gas employed in existing SAMs. 
Note, however, that gas particles are still converted to dissipation- 
less 'star' particles as dictated by the SAM. While cooling is rel- 
atively unimportant for the majority of the ICM, the low central 
entropy found in CC clusters cannot be reproduced in their simula- 
tions, so the large scatter towards the upper-luminosity edge of the 
observed local Lx-T relation is not recovered. This will be less of 
an issue at high-redshift since only a small fraction of clusters have 
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Vikhlinin etalll 20071) . due to the higher 
Jeltema et al]|2005h . 



cool cores at z > 0.5 (e.g. 
rate of major mergers (e.g 

A more self-consistent approach would be to fully couple a 
SAM to a radiative simulation, so that the gas distribution in the 
simulation is used to inform the SAM. This extension of the semi- 
analytic technique would require the simulation and the SAM to be 
coupled in such a way that both can be undertaken simultaneously. 
Extensive testing would be necessary to ensure that such a model 
was as successful as current S AMs in reproducing observed galaxy 
properties. Such a scheme is a long-term goal of our work but is 
beyond the scope of this paper. 

In order to elucidate the effect of feedback from galaxy forma- 
tion on the evolution of scaling laws, we compare the predictions of 
our feedback model with those of two other models in the Millen- 
nium Gas series. The first of these simulations incorporates gravi- 
tational heating only and is thus referred to as the GO run. This is 
useful as a base model. Given that the only source of gas entropy 
changes in the GO run is gravity, we would expect a self-similar 
cluster population to be form ed. This is generally found to be the 
case in such simulations (e .g. Navarr o et al ]|l995r,lEkeetalJll998l ; 



IVoit et al.ll2005l : lAscasibar et all2006l : IStanek et alj|2010h 

The second reference simulation includes high-redshift pre- 
heating and radiative cooling, in addition to shock heating. We 
name this the PC run. Preheating raises the entropy of the ICM be- 
fore gravitational collapse, preventing gas from reaching high den- 
sities in central cluster regions and thus reducing its X-ray emis- 
sivity. This effect is greater in lower-mass systems, breaking the 



bles observations dBialek et al.l 2001; 


Brighenti & Mathews[|200l|; 


Muanwong et al. 2002; Borgani et al. 


20021: iTornatore et alJbo03l: 


Boraani et al. 2005; Hartlev et al. 2008; Stanek et alj|2010l). How- 



ever, preheating fails to account for the observed scatter about the 
mean relations, particularly on group scales, and gene rates overly- 
large isentropic cores in low-mass systems (e.g. IPonman et al.l 
l2003l:|Prattetai]|2006h . 

The cosmological model adopted in all three Millennium 
Gas simulations is a spatially-flat ACDM model with parameters 
ft m ,o = 0.25, Q h , = 0.045, «a, = 0.75, h = 0.73, n s = 1 and 
o"8,o = 0.9. Here fi m ,o, fib.o an d £2a,o are the total matter, baryon 
and dark energy density parameters, respectively, h is the Hubble 
parameter Ho in units of 100 km s _1 Mpc -1 , n s is the spectral in- 
dex of primordial density perturbations, and cts.o is the rms linear 
density fluctuation within a sphere of comoving radius 8h~ L Mpc. 
The subscript signifies the value of a quantity at the present day. 
These cosmological parameters are the same as those used in the 
original Millennium simulation and are consistent with a combined 
analysis of the first-year W ilkinson Microwave Anisotropy Probe 
(WMAP) data dSpergel et afll2003l) and data from the Two-degree- 
Field Galaxy Redshift Survey dColless et al]|200ll) . However, there 
is some tension between the chosen parameter values, particularly 
n B and o~8,o, and tho se derived from the seven-year WMAP data 
dKomatsu etal1l2010h . 

We now describe further details of our simulati ons. The GO 
and PC runs have alrea dy been discussed elsewhere ( Hartl ev" st al.l 
l2008l : IStaneketal.ll2010l) . so here we only briefly summarise their 
properties, focusing our attention mainly on the new FO run. 



3.1 The GO and PC simulations 

Initial conditions for the GO and PC runs were created at a redshift 
Zi = 49 by displacing particles from a glass-like distribution, so as 
to form a random realisation of a density field with a ACDM linear 



powe r spectrum obtained from CMBFAST dSeliak & Zaldarriaga 
1996). The amplitudes and phases of the initial perturbations were 
chosen to match those of the Millennium simulation. In both 
cases, the simulation volume is a comoving cube of side length 
L = 500/i _1 Mpc, as in the Millennium Simulation, containing 
5 x 10 s dark matter particles of mass ttidm = 1.42xlO lo /i _1 M0, 
and 5 x 10 s gas particles of mass m gas = 3.12 x 10 9 /i _1 M Q . 
These are some of the largest hydrodynamical A-body simulations 
ever carried out. Although the mass resolution is approximately 20 
time coarser than the Millennium Simulation, over 95% of haloes 
formed with a mass greater than 10 13 /i _1 Mq are within 100/i -1 
kpc and 5% of their original position and mass, respectively. 

The massive ly parallel TreePM A-body/SPH code GADGET- 
2 dSpringelhOol) was used to evolve both sets of initial conditions 
to z = 0, with full particle data stored at 160 output redshifts. 
The Plummer-equivalent gravitational softening length was fixed 
at e = 100ft -1 kpc in comoving coordinates until z = 3, then 
fixed in physical coordinates thereafter. The softening is thus ap- 
proximately 4% of the mean interparticle spacing, which has been 
shown to be the optimal choice for hydrodynam ical simulations 
dThomas & Couchrnanlll992l : lBorgani et al.ll2006h - 

The simple model of preheating employ ed in the PC simula- 
tion is similar to that of iBorgani et alj (2002). Briefly, the entropy 
of every particle is raised to 200 keV cm 2 at z = 4, thus creating 
an entropy 'floor'. In addition to preheatin g, there is also radia- 
tive co oling based on the cooling function of Sutherland & Dopita 
dl993l) . assuming a fixed metallicity of 0.3 Z© (a good approxi- 
mation to the me an metallicity of the ICM out to at least z = 1; 
iTozzi et alj2003l) . Once the temperature of a gas particle drops be- 
low 2 x 10 4 K, the hydrogen density exceeds pn = 4.2 x 10" 27 
g cm -3 and the density contrast is greater than 100, then it is con- 
verted to a collisionless star particle. However, the preheating is so 
extreme that star formation is effectively terminated at z = 4, so 
that less than 2% of the baryons are locked-up in stars at z = 0. 



3.2 The FO simulation 

We have adopted a different approach for the FO simulation than 
for the other two Millennium Gas runs. Rather than simulating 
the entire Millennium volume, we decided instead to resimulate 
a sample of several hundred galaxy groups and clusters extracted 
from the original Millennium simulation. In this way, it will be less 
time/resourc e consuming to develop and test future extensions of 
the model of lShort & Thomas! d2009l) . 

There are three distinct components of the FO run: a dark 
matter-only resimulation of each region containing a cluster from 
our sample, semi-analytic galaxy catalogues built on the halo 
merger trees of these resimulations, and hydrodynamical resimula- 
tions of the same regions to track the energy injection from model 
galaxies. We now discuss each stage of the modelling process in 
turn. 



3.2. 1 Dark matter cluster resimulations 

Our cluster sample consists of 337 objects identified in the z = 
output of the Millennium Simulation, spanning a broad mass range: 
1.7xl0 13 /i _1 M sC Mv« «S 2.9xl0 15 /i~ 1 M . The sample was 
constructed as follows. All clusters with A/ V ir > 5 x 10 14 /i _1 M Q 
were selected while, at lower masses, a fixed number of clusters 
were chosen at random per logarithmic interval in virial mass. For 
each cluster in the sample, we extract a spherical region from the 
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z — Millennium snapshot that is centred on the object in ques- 
tion and has a radius equal to twice the cluster virial radius. The 
particles contained within this sphere are traced back to z\ = 127 
(the initial redshift of the Millennium Simulation) to identify the 
Lagrangian regi on from which th e object formed. Multi-mass ini- 
tial conditions dTormen et al. 19 93) are then ma de by following a 
procedure similar to that of Spri ngel et 

aD {2008), 

but with one ma- 
jor difference: we do not increase the mass resolution in the La- 
grangian region of interest, but degrade the resolution exterior to 
this region instead. Particles in the 'high-resolution' region then 
have (approximately) the same mass as in the parent Millennium 
Simulation, ttidm = 8.61 x 10 s h~ x M©, but more distant regions 
are sampled with progressively more massive particles. Our res- 
imulations can thus be thought of as 'desimulations'. A glass-like 
configuration was used for the unperturbed particle distribution in 
the high-resolution regions of all our initial conditions, and the ini- 
tial particle displacements were imprinted using the Zeldovich ap- 
proximation. The amplitudes and phases of the initial perturbations 
again match those of the Millennium Simulation. 

We adopt the same mass resolution as the Millennium Sim- 
ulation in the high-resolution regions of our initial conditions for 
two reasons. First, the properties of our resimulated clusters should 
then agree well with those of their counterparts in the original Mil- 
lennium Simulation, except for small perturbations induced by dif- 
ferent representations of the tidal field. Second, we will be able to 
construct more comprehensive galaxy catalogues, and thus a more 
detailed model for feedback from galaxies, than if we had used the 
coarser mass resolution of the other Millennium Gas runs. 

The initial conditions for each cluster in our sample are then 
evolved to the present day with GADGET-2. Raw particle data 
is stored at the 64 output redshifts of the Millennium Simula- 
tion. At each output time, dark matter haloes are identified as viri- 
alised groups of high-resolution particles using a parallel friends- 
of-friends (FOF) algorithm. We adopt a stan dard FOF linkin g 
length of 20% of the mean particle separation JDavis et al. 1985), 
and only save groups that contain at least 20 particles, so that the 
minimum halo mass is 1.7 x 10 10 /i _1 Mq. Gravitationally-bound 
substructures orbiting within these FOF halo es are then found wit h 
a parallel version of the SUBFIND algorithm I Spr ingel et al.l200ll) . 
The resulting catalogues of groups and subgroups are written out 
alongside the snapshot data. 

In the Millennium Simulation, the gravitational force law was 
softened isotropically on a fixed comoving scale of e — 5/i _1 
kpc (Plummer-equivalent), corresponding to approximately 2% of 
the mean interparticle spacing. We employ a different softening 
scheme in our cluster resimulations. The Plummer equivalent grav- 
itational softening length is fixed to e = 9.2ft -1 kpc in physical 
coordinates below z = 3, and to e = 36.8/i~ 1 kpc in comov- 
ing coordinates at higher redshifts (for the high-resolution parti- 
cles). The z = softening length is thus 4% of the mean particle 
separation. Our choice of a larger softening scale allows us to re- 
solve more low-mass subhaloes, because two-body heating effects 
are less important, so we can construct more detailed semi-analytic 
galaxy catalogues. 

It is important that the high-resolution region remains free 
from contamination by more massive boundary particles during the 
course of a resimulation. Since we are interested in the evolution of 
cluster properties within rsoo, we have checked whether there are 
any boundary particles in this region for every cluster used in this 
study (see Section [3~3l below for a description of how we construct 
cluster samples from our resimulations). Only three objects were 
found to contain boundary particles interior to rsoo, of which two 



were significantly contaminated. We discarded these two objects 
from our cluster samples. 

3.2.2 Feedback from a galaxy formation model 

Dark matter halo merger trees are constructed in post-processing 
for each resimulated region using the stored subhalo catalogues. 
This is done by exploiting the fact that each halo will have a 
unique d escendant in a hierar chical scenario of structure forma- 
tion; see ISpringel et ail (2005) for details of the procedure. Us- 
ing these merger trees, we have generated galaxy catalogues for 
all our resimulations by ap plying the Munich L-Galaxies SAM of 
|Pe Lucia & Blaizotlj2007l) . Galaxy properties are saved at the same 
64 redshifts as the simulation data. We adopt ed the same set of 
model parameters as lDe Lucia & BlaizoJ J2007t> . For a full descrip- 
tion of the p hysical processes in corp orated in L-Galaxies, we re fer 
the reader to lCroton etal] d2006h and lDe Lucia & Blaizol d2007h . 

The information contained within the semi-analytic galaxy 
catalogues enables us to calculate energy feedback from galaxies 
into the ICM. We only give an outline the procedure here, since a 
complete account is given in Short & Thomas) J2009I) . 

For each cluster in our sample, we first identify all model 
galaxies that lie within a distance r v j r of the centre of the clus- 
ter halo at z = 0, then use the galaxy merger trees to find all their 
progenitors. We only consider feedback from this subset of galax- 
ies since this is sufficient to correctly determine the properties of 
the ICM within rsoo, the region we are interested in. For each of 
these galaxies, we use its merger tree to compute the change in 
stellar mass, AM,, and mass accreted by the central black hole, 
AMbh, over the time interval At between successive model out- 
puts. Knowledge of AM, enables us to incorporate star formation 
in our simulations as described in the following section. 

We compute the energy imparted to intracluster gas by Type 
II SNe using the L-Galaxies supernova feedback model. In this 
model, the total amount of energy released by SNe in a time At 
is proportional to AM, . However, some of this energy is assumed 
to be used up reheating cold gas in the galactic disk. Any en- 
ergy remaining after reheating is used to eject gas from the halo 
into the surroundi ng medium, heating the ICM (see equation 20 of 
ICroton et all2006h . 

The model of AGN feedback we employ is that o flBower et al.l 

(2008). In this scheme, the heat energy input into the ICM by AGN 
over a time period At is the minimum of 



ASheat = E,AM BH C 



and 



AE hc 



eSMBH 



AE E 



(5) 



(6) 



where c is the speed of light, Ai^Edd is the energy released by a 
black hole accreting at the Eddington rate in a time At, and e r is 
the efficiency with which matter can be converted to energy near the 
event horizon. Following convention, we set e r = 0.1, which is ap- 
propriate for radi atively efficient accretion onto a non-rapidly spin- 
ning black hole ( Sha kura & Sunvaevll 19731) . The parameter esmbh 
is related to the structure of the accretion disk itself. At low accre- 
tion rates, the accretion disk is expected to be geometrically thick 
and advection dominated, enabling efficient jet production and ef- 
fective radio mode feedback. As the accretion rate is increased, it 
is thought that the vertical height of the disk eventual ly collapses, 
so mu ch more of the accretion energy is radiated away. Bow er et"ai1 
(2008) assume that this change in disk structure occurs once the ac- 
cretion rate reaches Mbh = esMBH-MEdd, where esmbh = 0.02, 
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leading to an upper limit on the amount of energy available for 
heating intracluster gas (equation[6}. 



3.2.3 Hydrodynamical cluster re simulations 

To track the effect of energy feedback from galaxies on the thermo- 
dynamical properties of the ICM, we couple the L-Galaxies SAM 
to hydrodynamical resimulations of our clusters. We use the same 
multi-mass initial conditions for these resimulations as described 
above, but we add gas particles with zero gravitational mass. This 
ensures that the dark matter distribution remains undisturbed by the 
inclusion of baryons, so that the halo merger trees used to gener- 
ate the semi-analytic galaxy catalogues will be the same. Although 
baryons can influence the s tructure and merger histories of dark 
Saro et al. 2008; Romano-Diaz et al.l 

are ne- 



mueiic 
: et alj 
et alJ Ii 



2009; _ _ . _. . 

20ld ; lDuffvetalj|20ld) . such effects 



haloes ([St anek 
l2009l ; |Pedrosir 

glected in all modern SAMs, which are founded on merger trees 
derived from the dark matter distribution. We are forced to follow 
this route so that we can use SAM input into our simulations. This 
approximation is unlikely to have any significant impact on the re- 
sults presented in this paper. 

The number of gas particles added to each set of initial con- 
ditions is chosen in such a way that their true mass would be the 
same as in the other two Millennium Gas runs. Note that gas parti- 
cles are only added to the high-resolution region. We include gas at 
a lower resolution than the dark matter simply to ease the compu- 
tational cost of our simulations. The resolution we have adopted is 
sufficient to obtain numerically-converg ed estimates of bulk clu ster 
properties for systems with T > 2 keV dShort & Thomasll 20091) . 

Each set of initial conditions is evolved from z\ = 127 to 
z = with a version of GADGET-2 that has been modified to 
accommodate gas part icles with zero gravitational mass. As in 
IShort&Thomaslfe009h . cooling processes are neglected. The same 
softening scheme as used in the dark matter resimulations is applied 
to the gas particles since they do not influence the gravitational dy- 
namics. Whenever an SPH calculation is required, we assign gas 
particles their true mass, so that gas properties are computed cor- 
rectly. Gas particles are also given their true mass for simulation 
data dumps, with the mass of the dark matter particles accord- 
ingly reduced to (1 — /b)mDM = 7.06 x 10 s hT 1 Mq, where 
/b = fib,o/^m,o = 0.18 is the mean cosmic baryon fraction in 
our cosmological model. 

Once an output redshift is reached, temporary 'galaxy' parti- 
cles are introduced throughout the high-resolution region at posi- 
tions specified by the SAM galaxy catalogue. For each galaxy, we 
know the increase in stellar mass and energy released by SNe and 
AGN since the last output (see above). We use this information to 
fo rm stars and heat gas in the vicinity of model galaxies as detailed 
in lShort & Thomas! ( 120091) . Briefly, the A7V star = AM„/m gas gas 
particles nearest to each galaxy are converted into collisionless star 
particles, where a stochastic method is used to ensure AA^tar is an 
integer. Once star formation is complete, the heat energy available 
from SNe and AGN is distributed amongst all gas particles con- 
tained within a sphere of radius r2oo centred on the galaxy, in such a 
way that each particle receives an equal entropy boost. If the galaxy 
is not the central galaxy of a FOF group, then L-Galaxies approxi- 
mates V200 by the radius of a sphere enclosing a mass equal to the 
product of the number of particles in the host subhalo and the parti- 
cle mass. In this case, we distribute the available heat energy within 
this radius instead. The total energy contributed by such galaxies is 
about an order of magnitude less than for central galaxies anyway. 
Following the injection of entropy, the galaxy particles are removed 



and the simulation continues until the next output time, when the 
process is repeated. Note that increasing the frequency with which 
energy is injected into ICM has a negligible impact on our results, 
because the time interval between our chosen 64 model outputs is 
always less than the galaxy halo dynamical time. 

3.3 Constructing cluster catalogues and profiles 

Cluster catalogues are generated at several redshifts for the three 
Millenniu m Gas simulations usin g a procedure similar to that em- 
ployed by Muan wong et al. U2002I) , which we now briefly describe. 

The first step is to identify gravitationally-bound groups of 
dark matter particles with the FOF algorithm. This was done on 
the fly in the FO run, and we have group catalogues stored at all 
28 output redshifts between z = 1.5 and z — for each resim- 
ulated cluster. For the GO and PC runs, FOF groups were iden- 
tified in post-processing, setting the linking length to be 10% of 
the mean interparticle separation. Only groups with 500 particles 
or more were kept, corresponding to a minimum halo mass of 
7.10 x 10 12 /i -1 Mq. We produced group catalogues for these two 
simulations at seven different redshifts: z = 1.5, 1.25, 1, 0.75, 0.5, 
0.25, 0. 

The spherical overdensity method is then used to construct 
cluster catalogues. Briefly, a sphere is grown about the most 
gravitationally-bound dark matter particle of each FOF group un- 
til radii are found that enclose mean overdensities of A v i T (z), 
A = 200, A = 500, A = 1000 and A = 2500, relative to 
the critical density p CI (z). In cases where clusters overlap, we only 
keep the object with the largest mass within r2soo- We also discard 
clusters with fewer than 1000 particles at each overdensity, which 
corresponds to a minimum cluster mass of 8.61 x 10 n fo -1 Mq in 
the FO run, and 1.73 x 10 13 /i _1 M in the GO and PC runs. 

During the cluster identification process we compute a variety 
of cluster properties, averaged within each choice of scale radius. 
The relevant properties for this work are the total mass, gas mass, 
temperature and X-ray luminosity. The measure of temperature 
we ad opt is the spectroscopic-like temperature T s i l lMazzotta et al.l 
2004). In the Bremsstrahlung regime (T > 2 keV), this temper- 
ature estimator has been shown to provide the closest match to 
the actual spectroscopic temperature, T spGC , obtained by fitting X- 
ray spectra of simulated clusters with a single-temperature plasma 
model. The X-ray luminosity is approximated by the bolomet- 
ric emission-wei ghted luminosity, assuming the cooling function 
of ISutherland & Dopiti] Jl993h and a fixed metallicity of 0.3 Z Q . 
Three-dimensional gas density, spectroscopic-like temperature and 
entropy profiles are also computed for all our clusters by averaging 
particle properties within spherical shells, centred on the minimum 
of the dark matter potential. 

In this paper we focus our attention on the evolution of X-ray 
scaling relations for massive clusters only, since it is the most mas- 
sive objects that are observed at high redshift. Previous numerical 
studies that use a smaller simulation volume than the Millennium 
volume have been unable to resolve sufficiently large numbers of 
massive clusters to investigate this in detail. Limiting our scope in 
this way has two additional benefits. First, the cooling times are 
longer in the central regions of massive clusters than in groups, so 
the lack of self-consistent cooling in our feedback model is less of 
an issue. Second, the X-ray emission will be predominantly thermal 
Bremsstrahlung, so the spectroscopic-like temperature provides an 
accurate estimate of the spectral temperature. 

We construct samples of massive clusters for our study from 
the cluster catalogues as follows. The starting point is to discard all 
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Table 1. Number of clusters in each of the Millennium Gas simulations as 
a function of redshift, once the mass/temperature cut appropriate for each 
scaling relation has been made (see text). 



Relation 


1.5 


Redshift 
1 


0.5 





GO simulation 










Yx-M, Tgi-Af, L x -M 


25 


145 


549 


1109 




15 


107 


441 


946 


PC simulation 










Y x -M, T A -M, Lx-M 


14 


102 


410 


881 


Lx-T s \ 


13 


93 


376 


838 


FO simulation 










Yx-M, T Bl -M, Lx-M 


18 


75 


148 


187 


Lx-T s \ 


15 


67 


139 


186 




FO 
PC 
GO 

CC (Arnaud et al. 2009) 
NCC (Arnaud et al. 2009) 



clusters whose average density is too low for r25oo to be defined. 
We do this to remove any low-mass objects that may have erro- 
neous properties due to being subclumps falling into more massive 
neighbouring systems. For scaling laws involving the total cluster 
mass (Yx-M, T B \-M and Lx-M), we then impose a mass cut of 
M 50 o > lO 14 ^" 1 M Q at all redshifts of interest. For the L x -T s i 
relation, a cut is made in T B \ instead to ensure completeness in T s i; 
only clusters with a spectroscopic-like temperature greater than that 
corresponding to a mass of M500 = 10 14 /i _1 Mq on the mean 
T B \-M relation are kept. Table[T]lists the number of clusters in our 
final samples as a function of redshift for each of the Millennium 
Gas simulations. 



4 RESULTS AND DISCUSSION 

We now use the samples of massive clusters extracted from our 
three simulations to investigate differences in evolutionary be- 
haviour that arise from adopting a plausible model for feedback 
from SNe and AGN, rather than simple preheating or gravitational 
heating. We first present results at z = as they will form the basis 
for measuring the evolution of the thermal properties of the ICM 
with redshift. 

Throughout the remainder of this work, all uncertainties are 
quoted at the 68% confidence level. 



4.1 Cluster profiles at z — 

Radial cluster profiles are more sensitive to the precise manner in 
which non-gravitational cooling and heating processes are imple- 
mented in numerical simulations than X-ray scaling laws. There- 
fore, we start by examining whether our model for feedback from 
galaxies is able to explain the temperature and entropy profiles 
of observed low-r edshift clusters. The observational dataset we 
use is REXCESS dBohringer et alj|2007l) . a representative sample 
of 33 local (z < 0.2) clu sters drawn from the REFLEX cata- 
logue (Bohringer et al. 2004), all of which have been observed with 
XMM-Newt on. Temperature pro files for the REXCESS clusters are 
presen ted in lArnaud et al.l {2009 ), and entropy profiles in lPratt et al.l 
d2010L hereafter PAP10). We choose to compare with REXCESS 
for three reasons. First, REXCESS clusters were selected in lumi- 
nosity only, thus ensuring no morphological bias, in such a way 
as to sample the X-ray cluster luminosity function in an optimal 
manner. Second, distances were optimised in REXCESS so that 



Figure 1. Mean spectroscopic-like temperature profiles, with lcr scatter, 
obtained from the Millennium Gas simulations. The light and dark shaded 
regions enclose the mean profiles, plus lcr sca tter, of CC and NCC clusters 
in the REXCESS sample I Arnaud et al. 2009), respectively. Only clusters 
with a mass M500 ^ 10 14 /i — 1 Mq are considered. 



rsoo falls well within the XMM-Newton field-of-view, increasing 
the precision of measurements at large radii. Third, the same defi- 
nition of rsoo is used as in this work. 

To facilitate a fair comparison with our simulated data, we 
only consider observed clusters with a mass M500 > W^h^MQ. 
We have also resc a led th e observational data to account for the fact 
that lArnaud et al] d2009l) and PAP10 assumed a slightly different 
cosmological model in their analysis. 

In the following, it will prove useful to divide the REX- 
CESS sample int o CC and non-cool-core (NCC) systems. As in 
IPratt et al.l j2009l hereafter PCA09), clusters are classified as CC 
systems if they have a central gas density E(z)~ 2 n c (0) > 4.8 x 
10-V /2 cm- 3 . 



4.1.1 Temperature profiles 

In Figure Q] we display the average spectroscopic-like temperature 
profile of clusters in the FO simulation. For comparison, we also 
show average profiles obtained from the reference GO and PC sim- 
ulations and the observational data o f lArnaud et al. (2009). We dis- 
card profile data at radii less than the gravitational softening length 
and only plot the average profile if there are 10 or more clusters in 
a given radial bin. All profiles have been normalised to the char- 
acteristic halo temperature, T500, computed from the self-similar 
model: 

G urnn M 500 ,„ 

J500 = — -: , V) 

I Kb rsoo 

where /1 = 0.59 is the mean molecular weight for a fully ionised 
gas of primordial composition and inn is the mass of a hydrogen 
atom. In the case of the observed profiles, T500 is calculated at the 
redshift of each individual cluster. With this scaling we would ex- 
pect cluster profiles to coincide in the pure gravitational heating 
scenario. 

Clusters formed in the GO simulation are clearly cooler 
than observed. This is because the spectroscopic-like tempera- 
ture estimate is biased low by the cool, low-entropy cores of 
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accreted subhaloes that a re prevalent in GO cluster haloes (e.g. 
iMathiesen & Evrar3l200ll) . The temperature profiles of individual 
clusters are flat in core regions, but we see a slight decline in the av- 
erage scaled profile at small radii, r < O.lrsoo- This is because we 
are only averaging over the profiles of the most massive clusters at 
such radii, since the gravitational softening length is a smaller frac- 
tion of rsoo for these objects, and the normalisation of the scaled 
temperature profiles decreases with increasing mass. In the self- 
similar model the normalisation of the scaled profiles should not 
depend on mass; the reason for the mass-dependence in our GO 
simulation is that more massive clusters are less concentrated than 
their low-mass counterparts, i.e. even the dark matter is not truly 
self-similar. 

In the PC and FO runs the source of non-gravitational heating 
is completely different, but the net effect of the entropy injection is 
the same: cool subclumps are erased and the average temperature 
of the intracluster gas increases. Both simulations lead to temper- 
ature profiles that provide a good overall match to the observed 
profiles of NCC clusters, being nearly isothermal in core regions 
r < 0.15r5oo- Again, the down-turn in the average scaled profiles 
visible in core regions arises because, for small values of r/r 500, 
we are taking an average of the scaled temperature profiles of the 
most massive objects only, which have a lower normalisation than 
those of less massive systems. 

Recall that we have chosen to neglect gas cooling processes in 
our FO simulations, so it is not surprising that we do not reproduce 
the gradual decline in temperature seen in the central parts of CC 
clusters, where the cooling time is short. Although the PC run does 
incorporate radiative cooling, it has a negligible effect on the gas 
temperature since the entropy of the ICM is raised before structure 
formation commences, preventing gas from reaching high densities 
in cluster cores and cooling efficiently. 

Fully self-consistent simulations with radiative cooling and 
stellar feedback typically predict temperature profiles with a sharp 
spike at small cluster-centric radii, followed by a rapid drop in 
temperature moving further into the c ore fe.g. lBorganieTaT]|2004 
iNagai et alj|2007l ; ISiiacki et alj|2007l) . This is due to the adiabatic 
compression of gas flowing in from cluster outskirts to compen- 
sate for the lack of pressure support caused by too much gas cool- 
ing out of the hot phase. Temperature profiles of this form clearly 
conflict with the smoo thly declining (flat) p r ofiles of observed CC 
(NCC) clusters (e. g. ISanderson et"aUl2006l ; | Vikhli nin et al. I l2006l : 
lArnaud et af]|2009l) . 

Even if a feedback mechanism is able to reproduce observed 
temperature profiles, this does not guarantee that radiative cool- 
ing has been balanced. For ex ample, the stellar feedback scheme 
employed by iKav et alj J20071) is capable of producing tempera- 
ture profiles that are in reasonable agreement with observational 
data. However, the resulting stellar fraction within clusters was 
found to be ~ 25%, far in excess of the observed value of ~ 10% 
( iBalogh et al. 2001; Lin et al. 2003l; lBalogh et alj2008h . 



In recent work, iFabian et al. I d2010h have included a sub-grid 
model for AGN feedback in hydrodynamical simulations. They 
found that the additional heating from AGN was insufficient to pre- 
vent overcooling in massive clusters, again leading to too high a star 
formation efficiency and sharply peaked temperature profiles. How- 
ever, on the scale of galaxy groups, their AGN feedback scheme 
was abl e to successfully regula te the thermal structure of the ICM; 
sec also lMcCarthv et alj < l2010l) . 



4.1.2 Entropy profiles 

The entropy of intracluster gas increases when heat energy is in- 
troduced, and decreases when radiative cooling carries heat energy 
away. Entropy profiles thus preserve a record of the ph ysical pro- 
cesses respo nsible for similarity breaking in clusters (e.g. lVoit et alj 
l2002ll2003l) . 

If shock heating were the only mechanism acting to raise the 
entropy of the gas, then analytical models based on spherical col- 
lapse predict that entropy sc ales with radius as K oc r 1 ' 1 out- 
side of central cluster regions jTozzi & NormanlboOll) . Cosmolog- 
ical simulations that only include gravitational heating give rise to 
slightly steeper entropy profiles in cluster outskirts: K oc r 1 ' 2 (e.g. 
IVoitetal.ll200l;lNagaietai]|2007l) . 

Observed profiles are also typically found to scale as K oc 
r ' at large cluster-centric radii, flattening in central regions 



(e.g. Ponman et al. 2003; Sun ct al] 120091 ; ICavagnolo etail [2009 ; 
Sand erson et alj 120091 PAP10). However, the precise radius at 
which this flattening occurs varies considerably, depending on such 
factors as the temperature (mass) of the system and whether it has 
a CC or a NCC. In particular, h otter, more massive o bjects have 
a higher mean core entropy (e.g. ICavagnolo et alj|2009h . and the 
profiles of NCC clusters flatten off at significant ly larger radii than 
those of CC clusters (e.g. lSanderson et al.ll2009L PAP10). 

Figure [2] compares the entropy profiles obtained from each 
of our three simulations with observational data from REXCESS 
(PAP10). For illustrative purposes, we also show the power-laws 
K oc r 11 and K oc r 1 ' 2 , assuming an arbitrary normalisation in 
both cases. We have scaled all entropy profiles by the 'virial' en- 
tropy, -K500 , defined as 



K500 = 



^o, 500 



(8) 



where n c ,5oo is the average electron density within rsoo, given by 
500/ b £(z) 2 pcr,o 



n e,500 — 



(9) 



and /i c = 1.14 is the mean molecular weight per free electron. Note 
that K500 depends only on the total halo mass, so is independent of 
the thermodynamic state of the gas. 

The entropy profiles of clusters extracted from the GO run 
are indeed well described by the power-law K oc r 1 ' 2 for r > 
0.15r5oo, agreeing with the results of previous studies. There is also 
very little scatter about the mean, indicating self-similar scaling. 
For radii interior to 0.15r5oo there is more diversity; some of the 
clusters have nearly isentropic cores while others show no signs of 
flattening. Compared to the observed entropy profiles of CC clus- 
ters, the profiles of objects in the GO run have a similar slope, at 
least for r > 0.15rsoo, but the normalisation is systematically too 
low. In the case of NCC clusters, it is evident that the GO model 
cannot explain the shallow profiles characteristic of these systems. 

Clusters formed in the PC and FO simulations have entropy 
profiles that are broadly consistent with the theoretical scaling 
K oc r 1 ' 1 at large radii r > rsoo- This supports the idea that 
gravity dominates the ICM thermodynamics in the outer regions 
of clusters. As we move in towards the core from 7-500, the slope 
decreases and the profiles flatten off, providing a fair match to the 
observed entropy profiles of NCC clusters. However, on average, 
both the PC and FO models predict entropy profiles with a shal- 
lower slope than those of NCC clusters in central regions, resulting 
in an overestimate of the core entropy. Note that the FO run yields 
entropy profiles that are slightly closer to the observed NCC cluster 
profiles than the PC run. 
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Figure 2. Mean entropy profiles, with lcr scatter, obtained from the Mil- 
lennium Gas simulations. The light and dark shaded regions enclose the 
mean profiles, plus lcr scatter, of observed CC and NCC clusters from 
REXCESS (PAP10), respectively. We only consider clusters with a mass 
M s00 ^W li h- 1 M e . 



Figure 3. Scaled entropy at rsoo as a function of total mass within rsoo for 
clusters in the Millennium Gas simulations. For clarity, we do not display 
individual data points for the GO ran, but instead show the best-fit relation 
and the typical dispersion about this relation. Observational data for CC and 
NCC clusters in the REXCESS sample (PAP10) are also shown for compar- 
ison, along with lcr error bars. Only clusters with M500 If 10 14 /i — 1 Mq 
are considered. 



In both the PC and FO runs, we see a drop in the average 
scaled entropy profiles at small values of r/rsoo, where we are av- 
eraging over the profiles of just the most massive systems. This 
is because the normalisation of cluster entropy profiles decreases 
with mass in these simulations. To demonstrate this, in Figure [3] 
we plot the scaled entropy at rsoo as a function of M500 for each 
of our simulations, along with observational data from REXCESS 
(PAP 10). For clarity, we do not plot individual points for the GO 
simulation, but instead show the best-fit relation in log-log space, 
and the typical dispersion about this relation. The gradient of the 
best-fit line is very close to zero, so the normalisation of the scaled 
entropy profiles is independent of mass for GO clusters, as in the 
self-similar model. By contrast, our PC and FO models predict that 
the scaled entropy at rsoo is a decreasing function of mass, imply- 
ing that non-gravitational heating affects the entropy structure of 
the ICM out to larger radii in lower-mass systems. This is consis- 
tent with the expectation that non-gravitational processes are more 
influential at the low-mass end of the cluster population. Massive 
clusters in the PC and FO runs have a similar entropy at rsoo to 
their GO counterparts. Note that the observational data points also 
appear to suggest that the scaled entropy at rsoo decreases with 
increasing mass, being scattered about the PC and FO model pre- 
dictions. 

The mass-dependence of the normalisation of scaled cluster 
entropy profiles in the PC and FO runs explains why we see a larger 
scatter about the average profile than in the GO run. The scatter 
about the mean profile in the PC run is similar to that found in ob- 
served profiles of NCC systems, but the FO run generates clusters 
with a wider range of entropy profiles, leading to a larger scatter 
than is observed. 

Neither the PC or the FO models are capable of reproducing 
the steeply declining entropy profiles seen in CC clusters. In the 
case of the FO run, this problem could potentially be overcome by 
including radiative cooling in our model, since cooling acts to lower 
the entropy in dense central regions where the gas cooling time is 



short. As we have said, cooling is included in the PC run, but it is 
curtailed at high redshift by the preheating. 

Self-consistent 7V-body/SPH simulations that incorporate 
cooling, star formation and associated feedback are able to produce 
entropy profiles that resemble those of CC clusters, with a normal- 
isation in the outer parts of clusters that is higher than predicted by 
pure gravitational heating, and a steep slope that remains roughly 
constant all the way into the core (e.g. iBorgani et alj|2002l 12004k 
iKav et al.l2004l;lKavl20"04l ; lKav et alj2007l) . However, this success 
is usually achieved at the expense of excessive star formation. Fur- 
thermore, such simulations fail to reproduce the observed entropy 
profiles of NCC systems. 

4.2 X-ray scaling relations at z — 

We now discuss whether our feedback model generates local X-ray 
scaling laws that are compatible with observations, focusing on the 
Yx-M, Tsi-M, Lx-M and Lx-T s i relations. In each case, we com- 
pare with the corresponding relation derived from the low-redshift 
REXCESS data by PCA09. Wherever possible, we explain any dif- 
ferences that arise using the knowledge gleaned from our discus- 
sion of cluster profiles. The data of PCA09 are particularly suitable 
for a comparison with our simulated cluster samples because they 
tabulate spectral temperatures and luminosities within rsoo (where 
rsoo is defined as in this work), and their luminosities are bolomet- 
ric. Note that we have rescaled the observational data to allow for 
the slightly different choice of cosmological parameters adopted by 
PCA09. 

For each set of cluster properties, (X, Y) = (M, Yx), 
(M, T s i), (M,Lx) and (T s i, Lx), we fit a power-law scaling re- 
lation of the form 

E{z) n Y = C (jX , (10) 
to our simulated data points by minimising \ 2 m l°g space. Here 
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Table 2. Best-fit parameters (with lcr errors) for the z = X-ray scaling 
relations obtained from the Millennium Gas simulations. 



100.0 f 



Relation 


Co 


a 


°"l°gio Y 


GO simulation 








Yx-M 


4.202 ± 0.070 


1.547 ±0.014 


0.087 


T Bl -M 


3.931 ± 0.057 


0.554 ±0.012 


0.076 


L x -M 


18.58 ± 0.53 


1.203 ± 0.024 


0.148 


Lx-lsi 


37.6 ± 1.5 


2.004 ± 0.040 


0.137 


PC simulation 








Yx-M 


5.622 ± 0.052 


1.7805 ±0.0079 


0.045 


T a i-M 


6.310 ± 0.031 


0.5512 ± 0.0042 


0.024 


L x -M 


5.549 ± 0.088 


1.842 ±0.013 


0.076 




4.563 ± 0.055 


3.297 ±0.020 


0.063 


FO simulation 








Yx-M 


5.757 ±0.069 


1.692 ±0.016 


0.048 


T sl -M 


6.333 ± 0.049 


0.521 ±0.010 


0.031 


Lx-M 


6.17 ±0.15 


1.777 ±0.033 


0.098 


Lx-Tsi 


4.99 ±0.12 


3.296 ±0.065 


0.104 



Co is the best-fitting normalisation of the relations, and a is the best-fitting 
slope; see equation jl0\ . o"i O g 10 Y is th e scatter about the mean relation as 
defined by equation jilt . 

X = 5 x W^h- 1 M if X = M and X = 6 keV if 
X = T sh The normalisation Co has units of 10 14 M keV, keV 
and 10 44 h- 2 erg s 1 for Y = Yx, T s \ and Lx, respectively. Best- 
fitting parameters a and Co for each relation are summarised in Ta- 
ble^ The factor E(z) n is included to remove the self-similar evo- 
lution predicted by equations 0-©, where the index n — —2/3, 
-2/3, -7/3 and -1 for the Yx-M, T sl -M, Lx-M and Lx-T sl 
relations, respectively. We include this scaling factor simply to 'ad- 
just' observational data to z = for comparison with our z = 
simulated clusters. This will only be a small effect for the redshift 
range (z < 0.2) probed by the REXCESS sample. 

Scatter in the relations, fi ogl0 y, is quanitfied via the rms de- 
viation of log 10 Y from the mean relation: 



1 N 









- log 10 Co 



N - 

(11) 

where N is the number of individual data points (X iy Yi). The scat- 
ter about each relation is also listed in Table[2] 



4.2.1 The Yx-M relation 

The Yx.-M relation is particularly important as both simulations 
and observations indicate that Yx is a low-scatter mass proxy, even 
in the presence of significant dynamical activity. Figure [4] shows 
the local Yx-M relation obtained from the the Millennium Gas FO 
simulation, along with the relations derived from the GO and PC 
simulations and the observational data of PCA09. We define Yx as 
the product of the gas mass inside rsoo and the spectroscopic-like 
temperature in the 0.15rsoo < r ^ rsoo region, for consistency 
with the observations. 

In the self-similar model we would expect the slope of the 
Yx-M relation to be a = 5/3 (equation^. The Yx-M relation ob- 
tained from the GO simulation is significantly shallower than both 
this and the observed relation. The slope of the PC Yx-M relation 
is consistent with the observed slope a ~ 1.8 from lArnaud et al] 
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Figure 4. Yx as a function of total mass within rsoo f° r z = clusters 
in the Millennium Gas simulations. We do not display data points from 
the GO run for clarity. Observational data for CC and NCC clusters from 
REXCESS (PCA09) is shown for comparison, along with lcr error bars. 



d2007h at the lcr level. On the other hand, the FO run yields a Yx-M 
relation that is shallower than observations suggest, with a slope 
closer to that expected from self- similar scaling. This is similar to 
the result a « 1.7 obtained bv lKravtsov et al J (12006) using simu- 
lations with radiative cooling, star formation and supernova feed- 
back. The predicted scatter about the mean relation is similar in the 
PC and FO runs, being about a factor of 2 less than in the GO sim- 
ulation. We note that Yx does not appear to be as tightly correlated 
with cluster mass as T s \ in any of our simulations. 

The fact that both the PC and FO Yx-M relations lie close to 
the self-similar prediction implies that Yx must be relatively unaf- 
fected by the non-gravitational heating in these models. In the case 
of the FO run, this is presumably because AGN feedback removes 
gas from the central regions of haloes, reducing the gas mass within 
rsoo, but this is offset by an increase in gas temperature caused by 
the continual injection of entropy from galaxy formation. Preheat- 
ing evacuates haloes and increases the temperature of the intraclus- 
ter gas at high redshift instead, but the eventual outcome is essen- 
tially the same. 

It is important to note that the masses of the REXCESS clus- 
ters were not calculated from gas temperature and density profiles 
using the equation of hydrostatic equilibrium. This is because, by 
construction, REXCESS contains objects with a wide variety of 
dynamical states, so hydrostatic equilibrium may be a poor approx- 
imation in some cases. PCA09 estim ate cluster masses using the 
Yx-M relation o f lArnaudet al] [2007 ) instead, so the observational 
data points shown in Figure|4]all lie on this r e lation by construction. 

The Yx-M relation of lArnaud et all {2007) was calibrated 
using clusters with hydrostatic mass estimates. Simulations have 
shown that the assumption of hydrosta tic equilibrium can bias such 
mass estimates low by ~ 10 — 20% dRasia et alj 20061: Kay et all 
2007| : [Nagai et alj|200l iBurns et alj|2008l : IPiffaretti & Valdarninil 
aljb 



200a iMeneghetti et alj|2009h . This is primarily because of addi- 
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Figure 5. Spectroscopic-like temperature as a function of total mass within 
r 500 f° r 2 = clusters in the Millennium Gas simulations. We also dis- 
play observational data for CC and NCC clusters in the REXCESS sample 
(PCA09). 



tional pressure support provided by subsonic bulk motions in the 
ICM. Further bias may be introduced if there is significant pres- 
sure support from no n-thermal compon ents such as cosmic rays 
and magnetic fields i lLagana et al1l2010h . Therefore, the masses of 
clusters in the REXCESS sample are also likely to be underesti- 
mated by the same amount. 

If the total mass is indeed biased low by ~ 10 — 20%, then 
)"500 W1 U t> e underestimated by ~ 3 — 7%. Consequently, the gas 
mass within r^oo will be biased low, because it is obtained by in- 
tegrating the density profile out to rsoo, and the temperature will 
be overestimated (Figure [TJ. The former effect will dominate the 
latter, so we also expect Y% to be underestimated. However, this 
is unlikely to fully compensate for the bias in the total mass esti- 
mate, so there may actually be a small offset between our PC and 
FO Yx-M relations and the observed relation once any hydrostatic 
mass bias is accounted for. 



4.2.2 The T A -M relation 

The T s \-M relations obtained from the GO, PC and FO runs at 
2 = are compared with observational data from REXCESS in 
Figure [5] 

For a given mass, clusters in the GO run are much cooler than 
observed, because the spectroscopic-like temperature estimate is 
dominated by the cool, dense cores of merging subhaloes. Differ- 
ences in the distribution of this substructure drive fluctuations in 
T a i, leading to a large scatter about the mean relation. The scatter 
is ~ 2.5 (3) times larger than in the FO (PC) run. The slope of 
the GO T s i-M relation is shallower than observed, and also signifi- 
cantly shallower than a — 2/3 expected fom the self-similar model 
(see equation [TJ. The reason for this is that concentration depends 
on cluster mass in our GO simulation: low-mass clusters are more 



concentrated, and thus hotter, than their high-mass counterparts, 
flattening the T B \-M relation relative to the self-similar prediction. 

Non-gravitational heating in the PC and FO runs raises the 
mean temperature of the ICM above that expected from gravita- 
tional heating alone (Figure [TJ. At a fixed mass, PC and FO clus- 
ters are thus hotter than their GO counterparts, with temperatures 
close to those of observed clusters. Both the PC and FO T s \-M 
relations have a similar normalisation and slope, even though the 
manner in which entropy is injected into the intracluster gas is 
completely differe nt in each c ase. T he slope a « 0.55 is close 
to that derived by iKav et af] J2007h from a simulation incorpo- 
rating gas cooling and stellar feedback. However, the REXCESS 
relation has a significantly steeper slope: a = 0.633 ± 0.032 or 
a = 0.622 ± 0.031 (G. W. Pratt, priv. comm.), depending on the 
fitting procedure adopted (see PCA09 for details). Note that the ob- 
served slope is consistent with the self-similar value, indicating that 
the T spcc -M relation is relatively insensitive to baryonic physics. 

Cooling processes are neglected in our FO model and ineffi- 
cient in the PC model, so systems with a CC are not formed in either 
simulation. If we only consider the NCC clusters in REXCESS, the 
resulting T spcc -M relation is shallower: a = 0.613 ± 0.022 or 
a = 0.617 ± 0.022 (G. W. Pratt, priv. comm.), which is slightly 
closer to, but still steeper than, the slope obtained from our PC and 
FO simulations. The intrinsic scatter (i.e. once measurement errors 
have been accounted for) about the REXCESS T spec -M relation 
for NCC clusters is cri ogl0 y = 0.025 ± 0.015 (G. W. Pratt, priv. 
comm.), regardless of the regression method used, which is con- 
sistent with the scatter about the mean relation obtained from our 
PC and FO simulations. However, we note that the observational 
scatter will be an underestimate of the true dispersion since the 
masse s of REXCESS clus ters were derived from the Yx-M rela- 
tion of lArnaud et alj 120071) assuming no intrinsic scatter about that 
relation. 

It appears that a better fit to the observed data could be ob- 
tained, at least for the NCC clusters, if the the observational data 
points in Figure [5] were shifted to the right by ~ 10% in mass. 
This is consistent with the level of bias expected from hydrostatic 
estimates of cluster mass. 

4.2.3 The Lx-M relation 

In Figure[6]we show the three local Millennium Gas Lx-M scaling 
relations, plus observational data from PCA09. 

The peak of the X-ray emission in GO clusters is unresolved 
in our simulation, so the computed luminosities are not trustwor- 
thy. We present the GO Lx-M relation merely to illustrate how 
dramatically the model fails in reproducing the observational data. 

The Lx-M relations obtained from the PC and FO simulations 
are both steeper than anticipated from pure gravitational heating. 
Both relations have a slope a « 1.8, whereas the slope measured 
by PCA09 is a = 1.81 ± 0.10 or a = 1.96 ± 0.11. In the case of 
the FO run, this departure from self-similar scaling arises because 
AGN feedback expels gas from central cluster regions, reducing 
the gas density and thus X-ray emissivity. This effect is stronger in 
less massive systems, steepening the Lx-M relation as observed. 
By contrast, similarity-breaking is accomplished in the PC model 
by boosting the entropy of the ICM before gravitational collapse 
commences, preventing gas from reaching high densities in cluster 
cores and lowering the X-ray luminosity. Nevertheless, both the PC 
and FO models yield almost identical Lx-M relations at z = 0. 

The large scatter towards the upper edge of the observed 
Lx-M relation is due to systems with a highly X-ray luminous CC. 
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Figure 6. Bolometric X-ray luminosity as a function of total mass within 
r 500 f° r z = clusters in the Millennium Gas simulations. Observational 
data for CC and NCC clusters from REXCESS (PCA09) is also shown. 

The PC and FO simulations cannot account for this scatter since 
neither model can reproduce the steeply declining entropy profiles 
of CC clusters (Figure©. Quantitatively, PCA09 measure the scat- 
ter about the mean Lx-M relation to be o"i ogl0 (y) = 0.166±0.026, 
which is about a factor of 2 greater than the dispersion about our 
mean PC and FO Lx-M relations, even without accounting for the 
intrinsic scatter about the Yx-M relation underpinning the REX- 
CESS cluster mass estimates. 

Removing the CC clusters from the REXCESS sample leads 
to a shallower slope, a = 1.705 ± 0.094 or a = 1.766 ± 0.093, 
and a reduction in the measured scatter, <r logl0 (y) = 0.094±0.017, 
both of which agree very well with the predictions of our FO model. 
By contrast the PC Lx-M relation appears to be slightly too steep 
with too little scatter. The difference in scatter predicted by our two 
non-gravitational heating models is attributable to the greater diver- 
sity in the entropy profiles of objects formed in the FO simulation 
(Figure©. 

There is an apparent offset between our PC and FO Lx-M re- 
lations and the observational data for NCC clusters. The magnitude 
of this offset is about ~ 10% in mass, which could be accounted for 
by bias in the observational mass estimates due to the assumption 
of hydrostatic equilibrium. 

4.2.4 The Lx-T s \ relation 

The local GO, PC and FO Lx-T B \ scaling relations are displayed 
in Figure© along with observational data from REXCESS. Again, 
the GO relation is only shown for illustrative purposes. 

The PC and FO simulations produce almost identical Lx-T a \ 
relations, with a steep slope a « 3.3 owing to the th e breaking 
of self -similarity induced by non-gravitational heating. I Kay et all 
(2007) found a similar slope using a simulation with a self- 
consistent stellar feedback scheme. For comparison, PCA09 find 
a = 2.70 ± 0.24 or a = 3.35 ± 0.32 for the REXCESS sample. 



Figure 7. Bolometric X-ray luminosity as a function of spectroscopic-like 
temperature for 2 = clusters in the Millennium Gas simulations. X- 
ray properties are calculated within rsoo- For comparative purposes, we 
plot observational data for CC and NCC clusters in the REXCESS sample 
(PCA09). 



The dispersion about the REXCESS Lx-T spcc relation is 
cr loglo(y) = 0.288 ± 0.050 or cr loglo(y) = 0.318 ± 0.059. This is 
roughly 3 (5) times larger than the scatter about the FO (PC) Lx-Tsi 
relation, because no X-ray cores are formed in either simulation. 
However, the PC and FO Lx-T s \ relations seem to provide a good 
fit to the NCC clusters in REXCESS. The observed relation for 
NCC clusters has a slope a = 2.89±0.21 or a = 3.06±0.19, with 
corresponding scatter a logl0 ( Y ) = 0.116 ± 0.025 or 0"i O g lo (Y") = 
0.124±0.030. This is consistent with the scatter about our FO rela- 
tion, but the PC relation is too tight, with a dispersion that is about 
of a factor of 2 less. 



4.3 Evolution of cluster profiles from z = 1.5 to z — 

We have demonstrated that our feedback model can reproduce the 
observed properties of massive low-redshift clusters reasonably 
well, apart from those with a CC. We have also seen that the z — 
properties of clusters formed in the feedback run can be replicated 
almost exactly with a simplistic preheating model, where the en- 
tropy of the ICM is raised impulsively at z = 4, rather than by 
continual heating from SNe and AGN. Consequently, we cannot 
discriminate between these two models using data from local ob- 
servations. 

We now investigate whether feedback from galaxy formation 
leads to significantly different evolutionary behaviour than simple 
preheating. In this way we may be able to break the low-redshift 
degeneracy of the two models. We begin by examining the evolu- 
tion of cluster profiles from z — 1.5 to z — 0; this will help us to 
understand the predicted evolution of X-ray scaling relations dis- 
cussed in Section |4~4l below. No attempt is made to compare our 
results with observations since the available data is, as yet, very 
limited at high redshift. 
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4. 3. 1 Results for the PC simulation 

In Figure[8]we show the evolution of gas density, spectroscopic-like 
temperature and entropy profiles for clusters in the PC run. In the 
top row, we plot the profiles of the 10 hottest systems (i.e. the most 
massive objects) at each redshift and, in the bottom row, profiles 
for clusters in a narrow temperature range 3 keV ^ T B \ ^ 4 keV. 
Again, we only keep profile data at radii greater than the gravita- 
tional softening length. We scale the density, temperature and en- 
tropy profiles by n e ,soo, T500 and K500, respectively; see equations 
<EJ — . With this scaling we would expect to see no evolution of 
cluster profiles in the self-similar model, where the ICM is only 
heated by gravitational processes. We have confirmed that this is 
indeed the case in our GO run. 

Focusing on the profiles of the hottest clusters, we see clear 
signs of evolution beyond the self-similar prediction, which can be 
understood as follows. Imposing a uniform entropy floor at z — 4 
boosts the entropy of the ICM more in core regions than at large 
radii. Gas is driven out from central cluster regions, flattening the 
density profile and increasing the normalisation of the temperature 
profile, relative to the prediction from gravitational heating alone. 
Note that the temperature must increase if the slope of the density 
profile decreases to maintain pressure support. 

After preheating, ejected gas is gradually reincorporated into 
descendant haloes as hierarchical growth proceeds. Since the gas 
density has already been lowered by the preheating, it does not de- 
crease as rapidly as in the gravitational heating scenario, so we see 
an increase relative to the average gas density n Ci 5oo as z — > 0. 
Likewise, the preheated gas has a higher temperature and entropy 
than if the only source of heating was gravity. Therefore, as gas 
is accreted back onto descendant haloes, compression and shock 
heating raise its temperature and entropy at a lesser rate than in 
the gravitational heating model. This explains why we see a drop 
in gas temperature and entropy with redshift relative to T500 and 
-fsTsoo, respectively. Given that there is no further non-gravitational 
heating of intracluster gas in the PC run, the high-redshift entropy 
injection will become increasingly 'diluted' with time, and cluster 
profiles will eventually resemble those obtained from a simulation 
with gravitational heating only. 

At any given redshift, several of the 10 hottest objects are 
likely to have undergone a recent major merger, which could poten- 
tially affect the shape, dispersion and evolution of cluster profiles. 
To investigate this, we first compute the substructure statistic 

S = I*"™"*'! , (12) 

?"500 

for each of the 10 most massive clusters at each redshift of inter- 
est. Here, x c is the location of the dark matter potential minimum, 
which we take to be the cluster centre, and x com is the centre of 
mass of the cluster gas, defined by 

^ i mi(xj — x c ) 

Xcom — X c , (13) 

where the sums are over all gas particles within rsoo- Systems un- 
dergoing a major merger will be dyna mically disturbed and will 
thus have a larger value of S. Following lKav et all l l2007l) . we say 
that a cluster is disturbed if S > 0.1, and relaxed otherwise. At 
each redshift, we have found that the shape and dispersion of the 
radial profiles shown in the top row of Figure [8]remain almost un- 
changed if we only consider relaxed clusters in the sample of the 
10 hottest systems. This signifies that our results are not affected 
by cluster mergers. 

The scaled profiles of clusters with 3 keV ^ T s \ ^ 4 keV 



evolve in similar way to those of the most massive objects, but their 
shape is different at each redshift. In particular, their density and en- 
tropy profiles are flatter. The reason for these differences is that we 
are now considering lower-mass clusters. Consequently, preheat- 
ing is more effective at removing gas from their shallower potential 
wells, modifying their thermodynamic properties out to larger radii. 



4.3.2 Results for the FO simulation 

The evolution of the scaled gas density, spectroscopic-like temper- 
ature and entropy profiles of clusters in the FO run is illustrated in 
Figure [9] 

The main point to note from the scaled profiles of the 10 
hottest clusters is that they do not evolve with redshift. This means 
that the gas density, temperature and entropy scale in the same 
way as predicted by the self-similar model: n B oc E{z) 2 , T si 
E(z) 2/3 M 2/3 and K oc M 2/3 / E(z) 2/S . Recall that the PC model 
predicts substantial evolution of scaled cluster profiles over the 
same redshift range, essentially because haloes are simply 'recov- 
ering' from the extreme preheating at z — 4. By contrast, energy 
feedback from galaxies is a continual process in the FO run, act- 
ing to reduce the gas density and increase the entropy in central 
regions. The fact that entropy is injected in such a way that the 
shape of the scaled cluster profiles does not change with redshift 
is presumably attributable to the self-regulatory nature of the feed- 
back loop in the galaxy formation model underpinning our simu- 
lation. This behaviour is not peculiar to the most massive objects; 
we find that the scaled profiles of lower-mass systems evolve in a 
self-similar fashion too. We have also checked that cluster mergers 
have a negligible effect on the shape, dispersion and evolution of 
radial profiles by following the procedure outlined in the previous 
section. 

Turning our attention to clusters with a temperature 3 keV ^ 
T B \ ^ 4 keV, we find that their scaled profiles do evolve, in the 
opposite sense to that predicted by the PC model. The explanation 
for this behaviour lies in the fact that we are considering a different 
set of clusters at each redshift. Consider the z = 1.5 progenitors of 
clusters with a temperature 3 keV ^ T s \ ^ 4 keV at z = 0. These 
objects are less massive than clusters with a temperature in the same 
range at z — 1.5. Accordingly, they will have been more affected 
by the non-gravitational heating from AGN, so the gas density will 
be lower and the entropy higher in central regions. Since the scaled 
profiles of individual clusters do not evolve in the FO run, these 
differences are preserved until z — 0, so we see apparent signs of 
evolution when comparing the profiles of clusters with 3 keV ^ 
T s i ^ 4 keV at z = with clusters of the same temperature at 
higher redshift. 

4.4 Evolution of X-ray scaling relations from z = 1.5 to 

z = 

We have seen that feedback from galaxy formation leads to dramat- 
ically different evolution of cluster profiles than high-redshift pre- 
heating. It follows that we should see differences in the evolution 
of X-ray scaling laws too. We now examine whether the evolution 
of the Yx-M, T s \-M, Ly^-M and Lx-7si relations predicted by our 
FO simulation is compatible with observational data, and whether 
the data prefers this model over simple preheating. The datasets 
we use are the low-reds hift REXCESS sample of PCA09, and the 
high-redshift sample o f iMaughan etafl fe008al hereafter MJF08). 
We choose to compare with the data of MJF08 for several reasons. 
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Figure 8. Evolution of gas density (first column), spectroscopic-like temperature (second column) and entropy (third column) profiles for clusters in the 
Millennium Gas PC run. The top row shows the evolution of the ten hottest systems, while the bottom row shows the evolution of clusters with a temperature 
in the range 3 keV < T A ^ 4 keV. 




Figure 9. Evolution of gas density, spectroscopic-like temperature and entropy profiles for clusters in the Millennium Gas FO run. The layout of the plots is 
identical to Figure[8] 



First, this dataset is one of the largest high-redshift X-ray-selected 
cluster samples currently available, consisting of 115 clusters ob- 
served with Chandra. Second, it covers a broader redshift range 
(0.1 < z < 1.3) than any other existing large sample. Third, tem- 
peratures and bolometric luminosities were derived in the aperture 
< r ^ rsoo, using the same definition of rsoo as in this work. 
This agrees with the way in which these properties were calculated 
for our simulated clusters. Nevertheless, we must be careful not to 
over-interpret any comparison of our simulated data with observa- 
tions since strong selection biases limit our ability to perform a sta- 



tistically meaningful comparison. We discuss this further in section 
|4~4~5l below. 

To account for the evolution of X-ray observables, we define 
scaling relations by 

E(z) n Y = C{z) (J^-J , (14) 

where all quantities are as in equation 1101 . except that the nor- 
malisation is now a function of redshift. For each relation, we fix 
the slope a to the z — value found previously (see Table [2}, 
and compute the normalisation at each redshift by minimising \ 2 
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Table 3. Best-fit parameters (with lcr errors) for the evolution of the nor- 
malisation of the X-ray scaling relations predicted by each Millennium Gas 
simulation. 



Relation 


Co 


P 


GO simulation 






y x -M 


4.222 ± 0.028 


-0.267 ±0.011 


Ta-M 


3.941 ± 0.031 


—0.335 ± 0.013 


Ly^-M 


19.53 ± 0.65 


—0.243 ± 0.055 




39.4 ± 1.4 


0.370 ± 0.058 


PC simulation 






y x -M 


5.96 ± 0.23 


-0.330 ±0.066 


T sl -M 


6.317 ±0.025 


0.0423 ± 0.0065 


L x -M 


6.18 ± 0.42 


-0.90 ±0.11 




5.34 ±0.52 


-1.77 ±0.16 


FO simulation 






Yx-M 


5.683 ± 0.038 


0.054 ± 0.014 


Ta-M 


6.180 ±0.042 


-0.249 ±0.014 


L x -M 


6.334 ± 0.046 


0.748 ±0.015 


Lx-T B \ 


5.85 ± 0.14 


0.760 ± 0.050 



Co is the best-fitting normalisation, and (3 is the best-fitting slope. /3 char- 
acterises the evolution of the normalisation; see equations )14t and H5\ . 

in log space. The self-similar model predicts that the slope of each 
relation will be independent of redshift. In our simulations, we see 
small fluctuations in the slope with redshift, but there is no system- 
atic variation, justifying our assumption of a fixed slope. 
A power-law of the form 

C(z) = Co(l + zf, (15) 

is then fit to the normalisation data to determine the parameters 
Co and /? (note that this may cause Co to change slightly from the 
z — value given in Table O. Best-fitting parameters for each 
relation are listed in Table [3] Since we have included the E(z) n 
factor in equation J 1 4b . then we would expect the slope /3 to be 
zero if clusters do indeed evolve self-similarly. If /3 < or /3 > 0, 
then we say there is negative or positive evolution, respectively. 
Note that some authors do not scale out the expected self-similar 
behaviour first, so their definition of negative/positive evolution has 
a different meaning to ours. This is one reason why care must be 
taken when comparing the results of different studies. 

4.4.1 The Yx-M relation 

Figure [Tol illustrates how the normalisation of the Yx-M relation 
evolves in each of the Millennium Gas simulations. The observa- 
tional data of PCA09 and MJF08 shown in the figure was plotted 
as follows. For each cluster in the two datasets, we computed C(z) 
using equation J 14b . a ssuming a was fi xed to the slope of the local 
Yx-M relation of lArnaud et all J2007I) : a = 1.825 ± 0.090. Note 
that MJF08 adopt the same definition of Yx as used in both this 
work and REXCESS. 

The normalisation of the Yx-M relation evolves in a negative 
sense in the PC run, which can be understood as follows. Preheating 
drives significant amounts of gas beyond rsoo at high redshift, even 
in the most massive of haloes. Given some cluster at z = 0, its gas 
mass at high redshift will thus be less in the PC run than expected 
from the pure gravitational heating scenario. On the other hand, 
Figure [8] shows that the spectroscopic-like temperature increases 
with redshift relative to the self-similar prediction. Since Yx is the 



a 

1.00 0.79 0.63 0.50 0.40 




- H h PC 

X X GO 

□ PCA09 
°' 4 - O MJF08 



i i i 

0.0 0.1 0.2 0.3 0.4 

logic 0+ z ) 

Figure 10. Normalisation of the Yx-M scaling relation as a function of 
redshift for each of the Millennium Gas simulations. Low-redshift observa- 
tional data from REXCESS (PCA09) and the high-redshift data of MJF08 
is shown for comparison, lcr error bars are also plotted for the observational 
data. 

product of these two quantities, these opposing evolutionary trends 
will cancel each other out to some degree. However, the net effect is 
a more rapid decrease in Yx with redshift than anticipated from the 
self-similar model. Since the total mass of a preheated cluster will 
decrease with redshift at a rate akin to the self-similar prediction, 
this implies that any cluster on the PC Yx-M relation at z — will 
follow a steeper trajectory towards the bottom-left of the Yx-M 
plane than expected. The drop in normalisation with redshift will 
then be larger than in the self-similar model, so we see negative 
evolution. 

Unlike preheating at high redshift, our model for feedback 
from galaxies produces a cluster population whose properties scale 
in a self-similar manner (see the top row of Figure |9). We would 
therefore expect to see no evolution of the normalisation of the FO 
Yx-M relation, relative to the prediction of the self-similar model. 
From Figure [10] we see that this is almost the case; the value of /3 
obtained from fitting the normalisation data is close to zero. There 
is slight positive evolution which arises because the slope of the 
relation, a « 1.7, is marginally steeper than the self-similar value. 
This effect is explained in full in Appendix lAl We note that the sig- 
nificant negative evolution seen in the GO model can be explained 
in the same way. 

4.4.2 The T sl -M relation 

The normalisation of the T s \-M relations obtained from the GO, 
PC and FO runs is shown as a function of redshift in Figure [TJJ 
Also shown is the observational data of PCA09 and MJF08, where 
we have assumed a fixed slope equal to that of the local REXCESS 
T spcc -M relation: a = 0.633 ± 0.032 (G. W. Pratt, priv. comm.). 

The PC model predicts slightly positive evolution of the nor- 
malisation. If we take a particular object that lies on the T s \-Ad rela- 
tion at z = 0, then we know from Figure[8]that its temperature will 
decrease less rapidly with redshift than in the gravitational heating 
scenario. However, its total mass will decrease at a similar rate with 
redshift as it is dominated by the dark matter. Therefore, the cluster 
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Figure 11. Normalisation of the T s \-M scaling relation as a function of 
redshift for each of the Millennium Gas simulations. We also display obser- 
vational data from PCA09 and MJF08. 

will move towards the bottom-left of the T s \-M plane, following a 
shallower trajectory than expected from self-similar scaling argu- 
ments. It follows that, at some z > 0, the normalisation will be 
higher than expected, i.e. positive evolution. Although the degree 
of evolution is small, recall that the slope of the PC T s \-M relation 
is substantially shallower than the self-similar value (see TableO. 
In Appendix |A] we show how this can induce negative evolution, 
implying that the positive evolution predicted by the PC model is 
actually stronger than it appears. 

The normalisation of the T s \-M relation obtained from the FO 
run evolves in the opposite sense, being significantly negative. A 
priori we would expect to find /3 ~ in this model since entropy is 
injected in such a way that the profiles of individual clusters scale 
self-similarly (Figure[9j. The cause of the apparent evolution is that 
the slope is considerably shallower than the self-similar prediction. 
Given this difference in slope, the argument presented in Appendix 
|A]implies that individual clusters must decrease in mass by about a 
factor of three between z = and z = 1 to account for the ~ 20% 
drop in normalisation seen in Figure[TT] We have checked that this 
is indeed the case in the FO simulation. 



4.4.3 The Lx-M relation 

In Figure [12] we show how the normalisation of the Lx-M rela- 
tion evolves in each of our three simulations. The observational 
data of PCA09 and MJF08 is also plotted, assuming a fixed slope 
a = 1.96 ±0.11, which is the slope of the low-redshift REXCESS 
Lx-M relation. 

The normalisation of the Lx-M relation evolves in a negative 
manner in the PC run. To explain this, consider some preheated 
cluster that lies on the z — relation. Recall from Figure[8]that the 
PC model predicts a drop in gas density, relative to the self-similar 
prediction, as redshift increases. Since the dominant contribution 
to the X-ray luminosity is the gas density in the Bremsstrahlung 
regime, then the luminosity at some higher redshift will be lower 
than predicted by the self-similar model. Given that the total mass 
of the cluster at this redshift will be close to the value expected 
from self-similar evolution, then we will see an apparent decrease 



Figure 12. Normalisation of the Lx-M scaling relation as a function of 
redshift for each of the Millennium Gas simulations. For comparative pur- 
poses, we plot low and high-redshift observational data from PCA09 and 
MJF08, respectively. 



in normalisation of the Lx-M relation relative to the self-similar 
prediction. This corresponds to negative evolution. We would ex- 
pect the evolution to appear stronger if the slope of the PC Lx-M 
relation matched the self-similar value, rather than being consider- 
ably steeper. 

The density and temperature profiles of clusters formed in the 
FO run evolve in a self-similar fashion, so the X-ray luminosity will 
also scale self-similarly. Since the growth rate of a cluster is gov- 
erned primarily by the dark matter dynamics, we would thus expect 
the normalisation of the Lx-M relation not to evolve once the pre- 
dicted self-similar behaviour has been factored out. However, this 
will only be the case if the slope of the relation matches the self- 
similar value. In reality, feedback from SNe and AGN establishes 
a steeper slope, a ~ 1.8. As discussed in Appendix lAl this depar- 
ture from self-similarity leads to apparent positive evolution. From 
Figure[T2] we see that the normalisation of the FO Lx-M relation 
does indeed evolve positively, with a ~ 40% increase in normali- 
sation relative to the self-similar model between z — and z = 1. 
This is consistent with the fact that the masses of clusters in the FO 
run decline by roughly a factor of three over this redshift range. 



4.4.4 The Lx-T Bl relation 

We show the normalisation of the three Millennium Gas Lx-T a i 
relations as a function of redshift in Figure[T3] To plot the observa- 
tional data of PCA09 and MJF08, we have fixed the slope to that of 
the local REXCESS relation: a = 3.35 ± 0.32. 

The PC run predicts negative evolution of the Lx-T s \ rela- 
tion as well, which is more pronounced than for the Lx-M rela- 
tion. This is because the temperature increases relative to the self- 
similar prediction with redshift (Figure [8}, whereas the total mass 
decreases at a similar rate. Over a given redshift interval, the nor- 
malisation of the Z/x-Tsi relation decreases more than that of the 
Lx-M relation, so we see a larger drop in normalisation relative 
to the self-similar model, implying stronger negative evolution. As 
before, the evolution of the PC Lx-T B \ relation will have been tem- 
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Figure 13. Normalisation of the Lx-T s \ scaling relation as a function of 
redshift for each of the Millennium Gas simulations. Data from the obser- 
vational studies of PCA09 and MJF08 is also shown. 



pered somewhat, because the slope is steeper than the self-similar 
value. 

On the other hand, our model for feedback from galaxies 
leads to an Lx-T B \ relation with a positively-evolving normalisa- 
tion. Given that the X-ray luminosity and spectroscopic-like tem- 
perature both scale self-similarly in this model, we should see no 
evolution relative to that expected from self-similar theory. How- 
ever, it is evident from Figure [T3l that, between z — and z = 1, 
the normalisation of the Lx-T B \ relation increases by ~ 50% com- 
pared to the self-similar prediction. Again, this evolution arises be- 
cause the slope of the Lx-T s \ relation obtained from the FO sim- 
ulation, a ~ 3.3, is much steeper than the self-similar prediction 
a = 2. With this difference in slope, the magnitude of the posi- 
tive evolution can be readily explained by following the argument 
outlined in Appendix [A] using the fact that the temperatures of in- 
dividual clusters in the FO run drop by about a factor of 2 between 
z = and z = 1. 

We note that iKav et alj J2007I) find negative evolution of the 
Lx-T B \ relation using a fully self-consistent simulation with radia- 
tive cooling, star formation and supernova feedback. Their work is 
directly comparable to ours (same choice of overdensity, same tem- 
perature definition, no core excision, etc.), so this indicates that in- 
cluding AGN feedback changes the way in which the Lx-T s i rela- 
tion evolves, possibly because of the different redshift-dependence 
of the two feedback mechanisms. 



4.4.5 Comparison with observations 

Energy feedback from galaxies leads to substantially different evo- 
lution of the Yx-M, Tsi-M, Lx-M and Lx-T s \ relations than uni- 
form preheating. We now discuss which of the PC and FO models, 
if either, is preferred by the high-redshift data of MJF08. 

We begin by noting that, as in REXCESS, the masses of clus- 
ters in the sample of MJF08 were estimated from a Yx-M relation. 
Since this relation was calibrated using clusters with hydrostatic 
mass estimates, the masses of their high-redshift clusters are also 
likely to be biased low by ~ 10 — 20%. Therefore, the observa- 



tional data points shown in Figures QT] and [12] should all be shifted 
down by ~ 5 — 10% and ~ 20 — 30%, respectively. 

Once we have applied this correction, we find that all four 
scaling relations obtained from the FO run evolve in a manner 
broadly consistent with the observational data at low to moderate 
redshifts, z < 0.5. In the case of the Lx-M and Lx-T s i relations, 
there are hints that the positive evolution predicted by our feedback 
model provides a better match to the data at these redshifts than the 
PC model, although the observed scatter is large. Both the PC and 
FO runs predict similar results for the Yx-M and T a \-M relations 
at z < 0.5, and it is not possible to distinguish between the two 
models with the observational data. 

At higher redshift, z > 0.5, the observational data for the 
Yx-M and Lx-M relations seems to follow an upward trend, 
consistent with the positive evolution expected from our feedback 
model. The data also suggests that the T spec -M relation evolves in 
a positive sense at these redshifts, but in this case the PC model 
provides a better description of the observed evolution than the 
FO model. The PC model also predicts negative evolution of the 
Lx-T s \ relation for z > 0.5, consistent with the observational data. 

It is clear from this discussion that it is difficult to deduce 
whether the data of MJF08 favours our feedback model over simple 
preheating. For each relation, it seems as if the observed normali- 
sation data cannot be well fit by a single power-law. For example, 
Figure [T3] suggests that the evolution of the Lx-T spec relation is 
approximately self-similar (or possibly slightly positive) up until 
z ~ 0.5, then negative thereafter. This could be a signature of a 
change in the evolutionary behaviour of clusters that is not repro- 
duced by any of our models or, more probably, it may simply be an 
artifact of selection effects instead. 

At low to moderate redshifts, most clusters in the heteroge- 
neous sample of MJF08 come from samples based on the ROSAT 
All-Sky Survey (RASS), which are wide and shallow. Their rela- 
tively high flux limit corresponds to an intermediate mass limit at 
low redshift, but a much higher mass limit at moderate redshift, 
thus falling on a steeper part of the mass function. Given the large 
scatter in the Lx-M relation, this means that the number of objects 
scattered from the left to the right of the mass limit will grow rela- 
tive to the number scattered in the other direction. This increasing 
bias towards luminous systems as we transition from low to mod- 
erate redshift may explain the 'hump' in the observational data at 
z ~ 0.3 apparent in Figures [12] and [T3] Unfortunately, it is hard 
to quantify this effect since the sample of MJF08 is not cleanly 
selected, so the selection function is unknown. 

At high redshift, we expect the data of MJF08 to be less af- 
fected by selection biases for two reasons. First, their high-redshift 
clusters come from narrow an d deep samples, suc h as the Wide 
Angle ROSAT Pointed Survey teurenin et ail |2007|) and the 400 
Square Degree ROSAT PSPC Survey dHorner et alJl2008D . whose 
lower flux limit corresponds to a high-redshift mass limit that falls 
on a flatter part of the mass function than the mass limit of RASS- 
based surveys at moderate redshift. This implies smaller bias given 
the same scatter in the Lx-M relation as at lower redshift. Second, 
there should be less scatter about the mean Lx-M relation at hig h 
redshift due to the absence of cool cores (e.g. lVikhlinin et al J2007h . 
which would further reduce any bias. However, given the remain- 
ing uncertainties on the selection biases and the limited number of 
high-redshift clusters in the sample of MJF08, we cannot draw firm 
conclusions about the nature of the evolution at z > 0.5. 

To summarise, it is fair to say that the quality of current X-ray 
data is insufficient to place robust constraints on theoretical mod- 
els. In addition to small numbers of high-redshift clusters and large 
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measurement errors, existing heterogeneous samples are plagued 
by strong selecti on biases which can imitate genuine evolution. As 
demonstrated bv lPacaud etaD d2007h . correctly modelling the full 
source-selection process is crucial for measuring the evolution of 
scaling laws. In the near future, the XCS will provide a large sample 
of X-ray-selected clusters (~ 500 objects) with < z < 1.5 that 
have been analysed in a consistent manner across the full redshift 
baseline. The survey selection function will be well monitored, al- 
lowing selection effects to be properly included when analysing 
the evolution of the X-ray scaling relations. With such a dataset it 
will hopefully become possible to discriminate between theoretical 
models such as our PC and FO models, providing us with a valuable 
insight into cluster astrophysics. 



5 SUMMARY AND CONCLUSIONS 

In this paper, we set out to investigate the evolution of galaxy clus- 
ter X-ray scaling relations using numerical simulations. The evo- 
lution of scaling laws is crucial for constraining cosmological pa- 
rameters with clusters surveys, and also offers a potentially pow- 
erful probe of the cooling and heating processes operating in clus- 
ters. Our main objective was to determine how including additional 
feedback from AGN in simulations affects the predicted evolution, 
and whether this is consistent with observations. Given that there 
is a substantial body of observational and theoretical evidence in- 
dicating that AGN are key in shaping the properties of galaxy clus- 
ters, it is clearly important to address this issue. However, all evo- 
lution studies to date have been based on simulations that only in- 
corporate feedback from star formation. 

The simulation we have used for our study - the FO run - 
is a new member of the Millennium Gas suite, presented for the 
first time here. The basic objective of the Millennium Gas Project 
is to add gas to the structures found in the original Millennium 
Simulation. The Millennium Gas simulations are ideal for study- 
ing the evolution of cluster properties, because their large volume 
(500 3 h~ 3 Mpc 3 ) means that we can resolve statistically significant 
numbers of massive clusters at all relevant redshifts. Furthermore, 
we can follow the formation of the richest clusters, which are the 
objects actually observed at high redshift. 

Feed back is implemented in our simulation using the hybrid 
scheme of lShort & Thomas! ( 2009), where the energy input into the 
ICM by SNe and AGN is calculated from a SAM of galaxy for- 
mation. This guarantees that feedback originates from a realistic 
galaxy population, whereas fully self-consistent simulations often 
predict excessive star formation on cluster scales. 

Our main conclusions can be summarised as follows. 

(i) Non-gravitational heating from SNe and AGN in the FO run 
produces a z — cluster population whose radial temperature and 
entropy profiles broadly agree with those of NCC clusters in the 
REXCESS sample (PAP10). In particular, the temperature profiles 
are close to isothermal in the core, and the entropy profiles are sig- 
nificantly flatter in central regions than the theoretical K oc r ' 
scaling observed in cluster outskirts. However, it seems that the en- 
tropy of the gas has been raised too much in the core, compared to 
the observational data. None of our clusters exhibit a gentle drop 
in temperature at small cluster-centric radii or a steadily declining 
entropy profile, both of which are characteristic of CC systems. 
This is because gas cannot lose entropy via radiative cooling in our 
simulation. We note that fully self-consistent hydrodynamical sim- 
ulations tend to suffer from the opposite problem, in the sense that 
radiative cooling leads to the over-production of cool cores. 



(ii) The Yx-M, T s i-M, Lx-M and Lx-T s \ scaling relations ob- 
tained from the FO run at z — generally match the local REX- 
CESS relations (PCA09), once we have accounted for the fact that 
the observed masses are likely to be biased low by ~ 10 — 20% due 
to the assumption of hydrostatic equilibrium. The exception is that 
we cannot explain the large scatter above the mean Lx-M and Lx- 
Tspoc relations seen in the observational data. This is because the 
source of this scatter is highly X-ray luminous CC systems which 
are not formed in our simulation. 

(iii) A crude model of non-gravitational heating from astrophys- 
ical sources in which the ICM is preheated at z — 4, rather than in 
response to galaxy formation, can produce a population of clusters 
whose z = properties closely resemble those of objects formed 
in the FO run. In fact, the two model predictions are so similar that 
they cannot be distinguished using high-quality local observations. 

(iv) Density, temperature and entropy profiles of individual clus- 
ters in the FO run all evolve in a self-similar fashion from z — 1.5 
to z — 0, although feedback from galaxies has modified their shape 
compared to that expected from pure gravitational heating. We sus- 
pect this is linked to the self-regulation of cooling and heating in 
the underlying model of galaxy formation. 

(v) The profiles of preheated clusters do not scale self-similarly. 
This is because the injection of entropy at high-redshift acts to re- 
move gas from central cluster regions, lowering the gas density 
and increasing its temperature. Following preheating, the proper- 
ties of the ICM can only be modified by gravitational processes, so 
the effect of the preheating is gradually erased and cluster profiles 
will eventually resemble those of clusters that have been subject to 
gravitational heating only. This 'recovery' from preheating is what 
drives the apparent evolution of cluster profiles relative to the self- 
similar model. 

(vi) Feedback from galaxy formation in our FO model leads 
to positive evolution of the Yx-M, Lx-M and Lx-T B \ relations, 
and negative evolution of the Tgi-M relation. By contrast, preheat- 
ing leads to sca ling relations that evolve in the opposite sense. 
iKav et al.l d2007t) also reported negative evolution of the Lx-T s i 
relation using a simulation with a self-consistent stellar feedback 
scheme. This suggests that additional heating from AGN feedback 
changes the way in which scaling laws evolve, possibly because 
AGN heating is still important in cluster cores at low-redshift, long 
after the peak of star formation. We have investigated whether the 
evolution predicted by our feedback model is consistent with X-ray 
observations of high-redshift clusters. Unfortunately, the large sam- 
ples of high-redshift clusters currently available are not cleanly se- 
lected, w hich is problematic since it may generate spurious evolu- 
tion (e.g. lPacaud et alj2007l) . This is possibly why different obser- 
vational studies give contradictory results. Consequently, we have 
not been able to decide whether our FO model provides a better de- 
scription of reality than simple preheating. However, it is encourag- 
ing that the evolutionary behaviour predicted by the two models is 
distinct, particularly in the case of the Lx-M and Lx-T s \ relations, 
so that we could potentiall y distinguish be tween them, and also 
other models (such as that o f lKav et al .12001) , when higher-quality 
data becomes available. As an example, the XCS will soon provide 
the largest ever sample of X-ray clusters selected with well-defined 
criteria, extending out to z ~ 1.5. Likewise, large high-redshift 
cluster samples are expected from SZ surveys currently underway. 
With such datasets, a rigorous comparison between theory and ob- 
servation will become possible, so that we will be able to use the 
evolution of cluster scaling laws as an additional constraint on mod- 
els of non-gravitational heating in clusters. 
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The Millennium Gas FO simulation introduced here is the 
only existing simulation that is large enough to follow the evolution 
of significant numbers of massive clusters at reasonable resolution, 
while also attempting to include some of the main physical pro- 
cesses involved in cluster formation: star formation and feedback 
from both SNe and AGN. Although our feedback model can gen- 
erally reproduce several key observational properties of clusters, at 
least for those without a cool core, it does have its limitations. In 
the fut ure, we plan to enhance the hybrid model of Sho rt & Thomas! 
(2009) in two major ways. 

First, we will adapt the model to follow the metal enrich- 
ment of intracluster gas by Type II and Type la SNe. This is im- 
portant since radial abundance profiles derived from X-ray ob- 
servations provide valuable constraints on the feedback mecha- 
nisms responsible for injecti n g metals into the di f fuse phase (e.g. 
De Grandi & Molendil 1200 ll ; iTamura et all |2004 IVikhlinin et all 
2005). Cora (2006J) ha ve already used a similar approach to tackle 
this problem (see also ICora et al. 2008). However, they neglected 
energy feedback from SNe and AGN in their hybrid model, which 
will have a significant impact on the way metals are distributed 
throughout the ICM. 

Second, we aim to self-consistently incorporate radiative cool- 
ing into the model as well, rather than relying on the simple cool- 
ing recipes employed in SAMs. These recipes usually assume that 
haloes have a spherically symmetric isothermal gas distribution but, 
in general, neither of these assumptions will hold in hydrodynam- 
ical simulations. To circumvent this problem, we intend to fully 
couple SAMs to radiative simulations, so that the gas distribution 
in the simulation governs star formation, black hole growth and 
associated feedback in the SAM. This is a non-trivial task, requir- 
ing the simulation and SAM to be run simultaneously. With this 
modification we hope to be able to reproduce the roughly bimodal 
distribution of core entropies found in real clusters. 
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Figure Al. Diagram to illustrate how evolution of a general X-ray scaling 
relation can be induced if the slope, a, differs from the self-similar value 
«ss- The arrow shows how the cluster highlighted by an open circle evolves 
with increasing redshift. See text for mathematical details. 
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APPENDIX A: MOCK EVOLUTION OF X-RAY SCALING 
RELATIONS 

In this Appendix we demonstrate how evolution of the scaling rela- 
tions can arise if the slope is different to the self-similar prediction, 
even if cluster properties scale self-similarly themselves, as in our 
FO simulation. We present our argument in terms of the evolution 
of a general scaling relation of the form d 14b . A diagram of the 
situation under consideration is shown in Figure lATI 

Consider a cluster at Z\ = that is located at the point (X\ , Y\) 
on a Y-X relation with a self-similar slope ass and normalisation 
Cf s , so 



(Al) 



Suppose the cluster evolves self-similarly until, at some redshift 
Zh > 0, it is located at the point (Xh, Yh), then 
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In reality, we expect non-gravitational cooling and heating 
processes to alter the slope of the Y-X relation from the self- 
similar prediction. Now suppose that our cluster was located at the 
same point, (Xi, Y\), at Z\, but lay on a Y-X relation with a slope 
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Evolution of cluster scaling relations 
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a where, without loss of generality, we take a < ass- If we as- 
sume self-similar scaling of cluster properties, as predicted by our 
feedback model, then the position of the cluster at Zh. will again be 
(Xh, ih). We can then write 

Y\ = Ci (^r\ , (A3) 



and 



A' 



where Ci and Ch are the normalisations of the Y-X relation at 
z\ and Zh, respectively. In the self-similar model we would have 
Ch/Ci = 1. From equations JA3b and dA4t it follows that 

a =E{zh) yAk) ' (A5) 



but equations dAlt and dA2t imply 

Yh 

y \Xi 

so 

n / v \ a ss- a 

Oh / -A_h 



= ) (A6) 
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Since Xh < A"i and ass > a, then Ch/C < 1, so we see a 
decrease in normalisation relative to the self-similar prediction with 
redshift, i.e. negative evolution. If the slope had been steeper than 
the self-similar value, a > ass, then we would have seen apparent 
positive evolution. 



